############################################################################################################################
########The wrong place at the wrong time? Territorial autonomy and conflict during regime transitions########
#####Global analyses (main and robustness checks)#####
#####Andreas Juon & Daniel Bochsler 2023 (Comparative Political Studies)#####
############################################################################################################################

####################################################################################
#########1. Libraries, data import, definitions#########
####################################################################################

#######1.1 Libraries#######
library(plyr)
library(dplyr)
library(zoo)
library(gridExtra)
library(grid)
library(glm.predict)
library(boot)
library(fastDummies)
library(stargazer)
library(car)
library(sensemakr)

#######1.2 Functions#######
get_confint<-function(model, vcovCL){
  t<-qt(.975, model$df.residual)
  ct<-coeftest(model, vcovCL)
  est<-cbind(ct[,1], ct[,1]-t*ct[,2], ct[,1]+t*ct[,2])
  colnames(est)<-c("Estimate","LowerCI","UpperCI")
  return(est)
}
cluster.se<-function(model, cluster)
{
  require(multiwayvcov)
  require(lmtest)
  vcovCL<-cluster.vcov(model, cluster)
  coef<-coeftest(model, vcovCL)
  return(coef)
}

#######1.3 Reading data#######
setwd("[!!!Set working directory containing the data file!!!]")
df_global <- read.csv("./juon_bochsler_cps_autonomy.csv",sep=",")
df_global <- subset(df_global, status_no < 3)#excluding groups with dominant or monopoly status
df_global$yearf <- as.factor(df_global$year)
df_global$region <- as.factor(df_global$region)
df_global <- dummy_cols(df_global, select_columns =c("region"))

#######1.4 Defining variable lists#######
group_controls <- c("state_control","included","size","log(distance_border+0.00001)","tek_sup_irre_any","tek_sdm","oil_area_pct")
country_controls <- c("log(gdppc_l1)","log(pop_l1)","sd_sum_l1","polityca")
sdm_years <- c("years_no_sd_l1", "I(years_no_sd_l1^2)", "I(years_no_sd_l1^3)")
violsdm_years <- c("sd_l1","years_no_violsd_l1", "I(years_no_violsd_l1^2)", "I(years_no_violsd_l1^3)")
additional_controls <- c("cleavage_mean","low_ratio","high_ratio")
vars_polityc <- c("autonomy","d5_tt_polityc","autonomy:d5_tt_polityc","d5_tt_polityc:autonomy","tt_polityc_shock","autonomy:tt_polityc_shock","tt_polityc_shock:autonomy")
vars_vdem2 <- c("autonomy","d5_tt_vdem2","autonomy:d5_tt_vdem2","d5_tt_vdem2:autonomy","tt_vdem2_shock","autonomy:tt_vdem2_shock","tt_vdem2_shock:autonomy")
vars_vdem2_shock1 <- c("autonomy","tt_vdem2_shock1","autonomy:tt_vdem2_shock1","tt_vdem2_shock1:autonomy")
vars_vdem2_shock5 <- c("autonomy","tt_vdem2_shock5","autonomy:tt_vdem2_shock5","tt_vdem2_shock5:autonomy")
vars_vdem2_type5 <- c("autonomy","d3_tt_vdem2","autonomy:d3_tt_vdem2","d3_tt_vdem2:autonomy","tt_vdem2","autonomy:tt_vdem2","tt_vdem2:autonomy")
var_vdem2_type3 <- c("autonomy1_inc0","autonomy1_inc1","autonomy0_inc1", "d5_tt_vdem2", "autonomy1_inc0:d5_tt_vdem2", "d5_tt_vdem2:autonomy1_inc0","autonomy1_inc1:d5_tt_vdem2","d5_tt_vdem2:autonomy1_inc1","autonomy0_inc1:d5_tt_vdem2","d5_tt_vdem2:autonomy0_inc1")
var_vdem2_type4 <- c("autonomy1_inc0","autonomy1_inc1","autonomy0_inc1", "tt_vdem2_shock:autonomy1_inc0","autonomy1_inc0:tt_vdem2_shock","tt_vdem2_shock2:autonomy1_inc0","autonomy1_inc0:tt_vdem2_shock2","tt_vdem2_shock:autonomy1_inc1","autonomy1_inc1:tt_vdem2_shock","tt_vdem2_shock2:autonomy1_inc1","autonomy1_inc1:tt_vdem2_shock2","tt_vdem2_shock:autonomy0_inc1","autonomy0_inc1:tt_vdem2_shock","tt_vdem2_shock2:autonomy0_inc1","autonomy0_inc1:tt_vdem2_shock2")
vars_vdem2_types <- c("d5_tt_vdem2","d5_tt_vdem2_st","d5_tt_vdem2_dem_or_aut","d5_tt_vdem2_dem","d5_tt_vdem2_aut","d3_tt_vdem2","d3_tt_vdem2_st","d3_tt_vdem2_dem_or_aut","d3_tt_vdem2_dem","d3_tt_vdem2_aut","tt_vdem2","tt_vdem2_st","tt_vdem2_dem_or_aut","tt_vdem2_dem","tt_vdem2_aut","autonomy")
group_controls_org <- c("state_control","included","size","log(distance_border+0.00001)","tek_sup_irre_any","tek_sdm","oil_area_pct","lorg_number_any_g")
country_controls_org <- c("log(gdppc_l1)","log(pop_l1)","sd_sum_l1","polityca")

#######1.5 Defining variable reporting order for stargazer#######
vars.order_main1 <- c("autonomy","sr_transition","autonomy:sr_transition","sr_transition:autonomy","sr_shock","sr_shock2","tt_polityc_shock","tt_polityc_shock2","autonomy:sr_shock","sr_shock:autonomy","autonomy:sr_shock2","sr_shock2:autonomy","autonomy:sr_shock2","sr_shock2:autonomy","autonomy_included0","autonomy_included0:sr_transition","sr_transition:autonomy_included0","autonomy_included0:sr_shock","autonomy_included0:sr_shock2","sr_shock:autonomy_included0","sr_shock2:autonomy_included0","autonomy_included1","autonomy0_included1","year_since89","year_since892","I(year_since892)","sd_onsets_sum_after1989","terraut_up_sum_after1989","sd_escalation_sum_after1989","sd_l1","sd_l1_cs","terraut_1989 * sd_l1","terraut_1989:sd_l1","sd_l1:terraut_1989","terraut_1989:sd_l1_cs","terraut_l1_below1936","sd_l1","terraut_1936 * sd_l1","terraut_1936:sd_l1","terraut_1936:sd_l1","terraut_l1 * sd_l1","terraut_l1:sd_l1","terraut_l1:sd_l1","recent_concession5","terraut_up_cs","terraut_1989 * sd_period_terraut_up_l1","terraut_1989:recent_concession5","sd_period_terraut_up_l1:terraut_1989","terraut_1989:terraut_up_cs","terraut_1936","terraut_1936:sd_period_terraut_up_l1","sd_l1:terraut_1936","terraut_l1","terraut_l1:sd_period_terraut_up_l1","sd_l1:terraut_l1","recent_concession","terraut_1936 * sd_period_terraut_up_l1","terraut_1936:sd_period_terraut_up_l1","sd_period_terraut_up_l1:terraut_1936")
vars.order_main2 <- c(gsub("log\\(|\\)","",group_controls),gsub("log\\(|\\)","",country_controls),additional_controls)


####################################################################################
#########2. Main models#########
####################################################################################

###a) models
sr1_vdem2_sd <- glm(as.formula(paste("sd_onset", paste(c("d5_tt_vdem2", "autonomy", "autonomy:d5_tt_vdem2","included","yearf","as.factor(region)",group_controls,country_controls, sdm_years), collapse = " + "), sep = " ~ ")), data=subset(df_global, sd_l1 == 0), family=binomial(link="logit"),na.action = na.exclude)
cse_sr1_vdem2_sd <- data.frame(cluster.se(sr1_vdem2_sd,as.factor(subset(df_global, sd_l1 == 0)$cowcode))[, 2])
sr1_vdem2_violsd <- glm(as.formula(paste("violsd_onset", paste(c("d5_tt_vdem2", "autonomy", "autonomy:d5_tt_vdem2","included","yearf","as.factor(region)",group_controls,country_controls, violsdm_years), collapse = " + "), sep = " ~ ")), data=subset(df_global, violsd_l1 == 0), family=binomial(link="logit"),na.action = na.exclude)
cse_sr1_vdem2_violsd <- data.frame(cluster.se(sr1_vdem2_violsd,as.factor(subset(df_global, violsd_l1 == 0)$cowcode))[, 2])
srg1_vdem2_sd <- glm(as.formula(paste("sd_onset", paste(c("tt_vdem2_shock*autonomy","included","yearf","as.factor(region)",group_controls,country_controls, sdm_years), collapse = " + "), sep = " ~ ")), data=subset(df_global, sd_l1 == 0), family=binomial(link="logit"),na.action = na.exclude)
cse_srg1_vdem2_sd <- data.frame(cluster.se(srg1_vdem2_sd,as.factor(subset(df_global, sd_l1 == 0)$cowcode))[, 2])
srg1_vdem2_violsd <- glm(as.formula(paste("violsd_onset", paste(c("tt_vdem2_shock*autonomy","included","yearf","as.factor(region)",group_controls,country_controls, violsdm_years), collapse = " + "), sep = " ~ ")), data=subset(df_global, violsd_l1 == 0), family=binomial(link="logit"),na.action = na.exclude)
cse_srg1_vdem2_violsd <- data.frame(cluster.se(srg1_vdem2_violsd,as.factor(subset(df_global, violsd_l1 == 0)$cowcode))[, 2])
stargazer(type="text",sr1_vdem2_sd, sr1_vdem2_violsd, srg1_vdem2_sd, srg1_vdem2_violsd,se=c(cse_sr1_vdem2_sd, cse_sr1_vdem2_violsd, cse_srg1_vdem2_sd, cse_srg1_vdem2_violsd), column.labels = c("Model 1", "Model 2", "Model 3", "Model 4"),omit=c("years_no","yearf","factor"),  dep.var.labels.include=TRUE, dep.var.labels = c("SDM onset","Violent SDM onset","SDM onset","Violent SDM onset"), style = "ajps", notes.append = TRUE, notes.align = "l", model.numbers = F, notes = "Country-clustered errors in parentheses. Cubic term for time dependence of dependent variable and year- and region-fixed effects included but not reported.", order=c(paste0("^", vars_vdem2, "$"), ":", paste0("^", vars.order_main2, "$")), covariate.labels = c("Autonomy","Transition (0-5)","Autonomy x transition (0-5)","Transition proximity","Autonomy x transition proximity","Most powerful","Included","Group size","TEK irredentism","TEK SDMs","Petroleum % in settlement area","Previous SDMs in country","Normalized polity score","Distance border (log)","GDP p.c. (log)","Population (log)","Ongoing SD"))
###b) Wald tests
sr1_vdem2_sd_jointsig <- linearHypothesis(sr1_vdem2_sd, "autonomy + d5_tt_vdem2:autonomy = 0", test = "Chisq", vcov = cluster.vcov(sr1_vdem2_sd, subset(df_global, sd_l1 == 0)$cowcode))
sr1_vdem2_sd_equality <- linearHypothesis(sr1_vdem2_sd, "autonomy = autonomy + d5_tt_vdem2:autonomy", test = "Chisq", vcov = cluster.vcov(sr1_vdem2_sd, subset(df_global, sd_l1 == 0)$cowcode))
sr1_vdem2_violsd_jointsig <- linearHypothesis(sr1_vdem2_violsd, "autonomy + d5_tt_vdem2:autonomy = 0", test = "Chisq", vcov = cluster.vcov(sr1_vdem2_violsd, subset(df_global, violsd_l1 == 0)$cowcode))
sr1_vdem2_violsd_equality <- linearHypothesis(sr1_vdem2_violsd, "autonomy = autonomy + d5_tt_vdem2:autonomy", test = "Chisq", vcov = cluster.vcov(sr1_vdem2_violsd, subset(df_global, violsd_l1 == 0)$cowcode))
srg1_vdem2_sd_jointsig <- linearHypothesis(srg1_vdem2_sd, "autonomy + tt_vdem2_shock:autonomy = 0", test = "Chisq", vcov = cluster.vcov(srg1_vdem2_sd, subset(df_global, sd_l1 == 0)$cowcode))
srg1_vdem2_sd_equality <- linearHypothesis(srg1_vdem2_sd, "autonomy = autonomy + tt_vdem2_shock:autonomy", test = "Chisq", vcov = cluster.vcov(srg1_vdem2_sd, subset(df_global, sd_l1 == 0)$cowcode))
srg1_vdem2_violsd_jointsig <- linearHypothesis(srg1_vdem2_violsd, "autonomy + tt_vdem2_shock:autonomy = 0", test = "Chisq", vcov = cluster.vcov(srg1_vdem2_violsd, subset(df_global, violsd_l1 == 0)$cowcode))
srg1_vdem2_violsd_equality <- linearHypothesis(srg1_vdem2_violsd, "autonomy = autonomy + tt_vdem2_shock:autonomy", test = "Chisq", vcov = cluster.vcov(srg1_vdem2_violsd, subset(df_global, violsd_l1 == 0)$cowcode))


####################################################################################
#########3. Other conditional factors (appendix B.1)#########
####################################################################################

##########3.1 Autonomy and government inclusion##########
###a) models
sr7a_vdem2_sd <- glm(as.formula(paste("sd_onset", paste(c("d5_tt_vdem2", "autonomy1_inc0","autonomy1_inc1","autonomy0_inc1", "autonomy1_inc0:d5_tt_vdem2","autonomy1_inc1:d5_tt_vdem2","autonomy0_inc1:d5_tt_vdem2","yearf","as.factor(region)","state_control","size","log(distance_border+0.00001)","tek_sup_irre_any","tek_sdm","oil_area_pct",country_controls, sdm_years), collapse = " + "), sep = " ~ ")), data=subset(df_global, sd_l1 == 0), family=binomial(link="logit"),na.action = na.exclude)
cse_sr7a_vdem2_sd <- data.frame(cluster.se(sr7a_vdem2_sd,as.factor(subset(df_global, sd_l1 == 0)$cowcode))[, 2])
sr7a_vdem2_violsd <- glm(as.formula(paste("violsd_onset", paste(c("d5_tt_vdem2", "autonomy1_inc0","autonomy1_inc1","autonomy0_inc1", "autonomy1_inc0:d5_tt_vdem2","autonomy1_inc1:d5_tt_vdem2","autonomy0_inc1:d5_tt_vdem2","yearf","as.factor(region)","state_control","size","log(distance_border+0.00001)","tek_sup_irre_any","tek_sdm","oil_area_pct",country_controls, violsdm_years), collapse = " + "), sep = " ~ ")), data=subset(df_global, violsd_l1 == 0), family=binomial(link="logit"),na.action = na.exclude)
cse_sr7a_vdem2_violsd <- data.frame(cluster.se(sr7a_vdem2_violsd,as.factor(subset(df_global, violsd_l1 == 0)$cowcode))[, 2])
srg7a_vdem2_sd <- glm(as.formula(paste("sd_onset", paste(c("tt_vdem2_shock*autonomy1_inc0","tt_vdem2_shock*autonomy1_inc1","tt_vdem2_shock*autonomy0_inc1","yearf","as.factor(region)","state_control","size","log(distance_border+0.00001)","tek_sup_irre_any","tek_sdm","oil_area_pct",country_controls, sdm_years), collapse = " + "), sep = " ~ ")), data=subset(df_global, sd_l1 == 0), family=binomial(link="logit"),na.action = na.exclude)
cse_srg7a_vdem2_sd <- data.frame(cluster.se(srg7a_vdem2_sd,as.factor(subset(df_global, sd_l1 == 0)$cowcode))[, 2])
srg7a_vdem2_violsd <- glm(as.formula(paste("violsd_onset", paste(c("tt_vdem2_shock*autonomy1_inc0","tt_vdem2_shock*autonomy1_inc1","tt_vdem2_shock*autonomy0_inc1","yearf","as.factor(region)","state_control","size","log(distance_border+0.00001)","tek_sup_irre_any","tek_sdm","oil_area_pct",country_controls, violsdm_years), collapse = " + "), sep = " ~ ")), data=subset(df_global, violsd_l1 == 0), family=binomial(link="logit"),na.action = na.exclude)
cse_srg7a_vdem2_violsd <- data.frame(cluster.se(srg7a_vdem2_violsd,as.factor(subset(df_global, violsd_l1 == 0)$cowcode))[, 2])
stargazer(type="text",sr7a_vdem2_sd, sr7a_vdem2_violsd, srg7a_vdem2_sd, srg7a_vdem2_violsd,se=c(cse_sr7a_vdem2_sd, cse_sr7a_vdem2_violsd, cse_srg7a_vdem2_sd, cse_srg7a_vdem2_violsd), column.labels = c("Model 1", "Model 2", "Model 3", "Model 4"),omit=c("years_no","yearf","factor"),  dep.var.labels.include=TRUE, dep.var.labels = c("SDM onset","Violent SDM onset","SDM onset","Violent SDM onset"), style = "ajps", notes.append = TRUE, notes.align = "l", model.numbers = F, notes = "Country-clustered errors in parentheses. Cubic term for time dependence of dependent variable and year- and region-fixed effects included but not reported.", order=c(paste0("^", var_vdem2_type3, "$"), ":", paste0("^", var_vdem2_type4, "$"), ":", paste0("^", vars.order_main2, "$")), covariate.labels = c("Autonomy, not included","Autonomy, included","Included without autonomy","Transition (0-5)","Autonomy, not included x transition (0-5)","Autonomy, included x transition (0-5)","Included without autonomy x transition (0-5)","Autonomy, not included x transition proximity","Autonomy, included x transition proximity","Included without autonomy x transition proximity","Most powerful","Group size","TEK irredentism","TEK SDMs","Petroleum % in settlement area","Previous SDMs in country","Normalized polity score","Transition proximity","Distance border (log)","GDP p.c. (log)","Population (log)","Ongoing SD"))
###b) Wald tests
sr7a_vdem2_sd_jointsig1 <- linearHypothesis(sr7a_vdem2_sd, "autonomy1_inc0 + d5_tt_vdem2:autonomy1_inc0 = 0", test = "Chisq", vcov = cluster.vcov(sr7a_vdem2_sd, subset(df_global, sd_l1 == 0)$cowcode))
sr7a_vdem2_sd_jointsig2 <- linearHypothesis(sr7a_vdem2_sd, "autonomy1_inc1 + d5_tt_vdem2:autonomy1_inc1 = 0", test = "Chisq", vcov = cluster.vcov(sr7a_vdem2_sd, subset(df_global, sd_l1 == 0)$cowcode))
sr7a_vdem2_sd_equality1 <- linearHypothesis(sr7a_vdem2_sd, "autonomy1_inc0 = autonomy1_inc0 + d5_tt_vdem2:autonomy1_inc0", test = "Chisq", vcov = cluster.vcov(sr7a_vdem2_sd, subset(df_global, sd_l1 == 0)$cowcode))
sr7a_vdem2_sd_equality2 <- linearHypothesis(sr7a_vdem2_sd, "autonomy1_inc1 = autonomy1_inc1 + d5_tt_vdem2:autonomy1_inc1", test = "Chisq", vcov = cluster.vcov(sr7a_vdem2_sd, subset(df_global, sd_l1 == 0)$cowcode))
sr7a_vdem2_violsd_jointsig1 <- linearHypothesis(sr7a_vdem2_violsd, "autonomy1_inc0 + d5_tt_vdem2:autonomy1_inc0 = 0", test = "Chisq", vcov = cluster.vcov(sr7a_vdem2_violsd, subset(df_global, violsd_l1 == 0)$cowcode))
sr7a_vdem2_violsd_jointsig2 <- linearHypothesis(sr7a_vdem2_violsd, "autonomy1_inc1 + d5_tt_vdem2:autonomy1_inc1 = 0", test = "Chisq", vcov = cluster.vcov(sr7a_vdem2_violsd, subset(df_global, violsd_l1 == 0)$cowcode))
sr7a_vdem2_violsd_equality1 <- linearHypothesis(sr7a_vdem2_violsd, "autonomy1_inc0 = autonomy1_inc0 + d5_tt_vdem2:autonomy1_inc0", test = "Chisq", vcov = cluster.vcov(sr7a_vdem2_violsd, subset(df_global, violsd_l1 == 0)$cowcode))
sr7a_vdem2_violsd_equality2 <- linearHypothesis(sr7a_vdem2_violsd, "autonomy1_inc1 = autonomy1_inc1 + d5_tt_vdem2:autonomy1_inc1", test = "Chisq", vcov = cluster.vcov(sr7a_vdem2_violsd, subset(df_global, violsd_l1 == 0)$cowcode))
srg7a_vdem2_sd_jointsig1 <- linearHypothesis(srg7a_vdem2_sd, "autonomy1_inc0 + tt_vdem2_shock:autonomy1_inc0 = 0", test = "Chisq", vcov = cluster.vcov(srg7a_vdem2_sd, subset(df_global, sd_l1 == 0)$cowcode))
srg7a_vdem2_sd_jointsig2 <- linearHypothesis(srg7a_vdem2_sd, "autonomy1_inc1 + tt_vdem2_shock:autonomy1_inc1 = 0", test = "Chisq", vcov = cluster.vcov(srg7a_vdem2_sd, subset(df_global, sd_l1 == 0)$cowcode))
srg7a_vdem2_sd_equality1 <- linearHypothesis(srg7a_vdem2_sd, "autonomy1_inc0 = autonomy1_inc0 + tt_vdem2_shock:autonomy1_inc0", test = "Chisq", vcov = cluster.vcov(srg7a_vdem2_sd, subset(df_global, sd_l1 == 0)$cowcode))
srg7a_vdem2_sd_equality2 <- linearHypothesis(srg7a_vdem2_sd, "autonomy1_inc1 = autonomy1_inc1 + tt_vdem2_shock:autonomy1_inc1", test = "Chisq", vcov = cluster.vcov(srg7a_vdem2_sd, subset(df_global, sd_l1 == 0)$cowcode))
srg7a_vdem2_violsd_jointsig1 <- linearHypothesis(srg7a_vdem2_violsd, "autonomy1_inc0 + tt_vdem2_shock:autonomy1_inc0 = 0", test = "Chisq", vcov = cluster.vcov(srg7a_vdem2_violsd, subset(df_global, violsd_l1 == 0)$cowcode))
srg7a_vdem2_violsd_jointsig2 <- linearHypothesis(srg7a_vdem2_violsd, "autonomy1_inc1 + tt_vdem2_shock:autonomy1_inc1 = 0", test = "Chisq", vcov = cluster.vcov(srg7a_vdem2_violsd, subset(df_global, violsd_l1 == 0)$cowcode))
srg7a_vdem2_violsd_equality1 <- linearHypothesis(srg7a_vdem2_violsd, "autonomy1_inc0 = autonomy1_inc0 + tt_vdem2_shock:autonomy1_inc0", test = "Chisq", vcov = cluster.vcov(srg7a_vdem2_violsd, subset(df_global, violsd_l1 == 0)$cowcode))
srg7a_vdem2_violsd_equality2 <- linearHypothesis(srg7a_vdem2_violsd, "autonomy1_inc1 = autonomy1_inc1 + tt_vdem2_shock:autonomy1_inc1", test = "Chisq", vcov = cluster.vcov(srg7a_vdem2_violsd, subset(df_global, violsd_l1 == 0)$cowcode))

##########3.2 Autonomy and prior SD violence##########
###a) models
sr7b_vdem2_sd <- glm(as.formula(paste("sd_onset", paste(c("d5_tt_vdem2", "autonomy1_prevv0","autonomy1_prevv1", "autonomy1_prevv0:d5_tt_vdem2","autonomy1_prevv1:d5_tt_vdem2","violsd_prev_l1","violsd_prev_l1:d5_tt_vdem2","yearf","as.factor(region)",group_controls,country_controls, sdm_years), collapse = " + "), sep = " ~ ")), data=subset(df_global, sd_l1 == 0), family=binomial(link="logit"),na.action = na.exclude)
cse_sr7b_vdem2_sd <- data.frame(cluster.se(sr7b_vdem2_sd,as.factor(subset(df_global, sd_l1 == 0)$cowcode))[, 2])
sr7b_vdem2_violsd <- glm(as.formula(paste("violsd_onset", paste(c("d5_tt_vdem2", "autonomy1_prevv0","autonomy1_prevv1", "autonomy1_prevv0:d5_tt_vdem2","autonomy1_prevv1:d5_tt_vdem2","violsd_prev_l1","violsd_prev_l1:d5_tt_vdem2","yearf","as.factor(region)",group_controls,country_controls, violsdm_years), collapse = " + "), sep = " ~ ")), data=subset(df_global, violsd_l1 == 0), family=binomial(link="logit"),na.action = na.exclude)
cse_sr7b_vdem2_violsd <- data.frame(cluster.se(sr7b_vdem2_violsd,as.factor(subset(df_global, violsd_l1 == 0)$cowcode))[, 2])
srg7b_vdem2_sd <- glm(as.formula(paste("sd_onset", paste(c("tt_vdem2_shock*autonomy1_prevv0","tt_vdem2_shock*autonomy1_prevv1","tt_vdem2_shock*violsd_prev_l1","yearf","as.factor(region)",group_controls,country_controls, sdm_years), collapse = " + "), sep = " ~ ")), data=subset(df_global, sd_l1 == 0), family=binomial(link="logit"),na.action = na.exclude)
cse_srg7b_vdem2_sd <- data.frame(cluster.se(srg7b_vdem2_sd,as.factor(subset(df_global, sd_l1 == 0)$cowcode))[, 2])
srg7b_vdem2_violsd <- glm(as.formula(paste("violsd_onset", paste(c("tt_vdem2_shock*autonomy1_prevv0","tt_vdem2_shock*autonomy1_prevv1","tt_vdem2_shock*violsd_prev_l1","yearf","as.factor(region)",group_controls,country_controls, violsdm_years), collapse = " + "), sep = " ~ ")), data=subset(df_global, violsd_l1 == 0), family=binomial(link="logit"),na.action = na.exclude)
cse_srg7b_vdem2_violsd <- data.frame(cluster.se(srg7b_vdem2_violsd,as.factor(subset(df_global, violsd_l1 == 0)$cowcode))[, 2])
stargazer(type="text",sr7b_vdem2_sd, sr7b_vdem2_violsd, srg7b_vdem2_sd, srg7b_vdem2_violsd,se=c(cse_sr7b_vdem2_sd, cse_sr7b_vdem2_violsd, cse_srg7b_vdem2_sd, cse_srg7b_vdem2_violsd), column.labels = c("Model 1", "Model 2", "Model 3", "Model 4"),omit=c("years_no","yearf","factor"),  dep.var.labels.include=TRUE, dep.var.labels = c("SDM onset","Violent SDM onset","SDM onset","Violent SDM onset"), style = "ajps", notes.append = TRUE, notes.align = "l", model.numbers = F, notes = "Country-clustered errors in parentheses. Cubic term for time dependence of dependent variable and year- and region-fixed effects included but not reported.", order=c(paste0("^", var_vdem2_type3, "$"), ":", paste0("^", var_vdem2_type4, "$"), ":", paste0("^", vars.order_main2, "$")), covariate.labels = c("Transition (0-5)","Autonomy, no previous violence x transition (0-5)","Autonomy, previous violence x transition (0-5)", "Previous violence x transition (0-5)","Autonomy, no previous violence x transition proximity","Autonomy, previous violence x transition proximity","Previous violence x transition proximity","Most powerful","Included","Group size","TEK irredentism","TEK SDMs","Petroleum % in settlement area","Previous SDMs in country","Normalized polity score","Transition proximity","Autonomy, no previous violence","Autonomy, previous violence","Previous violence","Distance border (log)","GDP p.c. (log)","Population (log)","Ongoing SD") )
###b) Wald tests
sr7b_vdem2_sd_jointsig1 <- linearHypothesis(sr7b_vdem2_sd, "autonomy1_prevv0 + d5_tt_vdem2:autonomy1_prevv0 = 0", test = "Chisq", vcov = cluster.vcov(sr7b_vdem2_sd, subset(df_global, sd_l1 == 0)$cowcode))
sr7b_vdem2_sd_jointsig2 <- linearHypothesis(sr7b_vdem2_sd, "autonomy1_prevv1 + d5_tt_vdem2:autonomy1_prevv1 = 0", test = "Chisq", vcov = cluster.vcov(sr7b_vdem2_sd, subset(df_global, sd_l1 == 0)$cowcode))
sr7b_vdem2_sd_equality1 <- linearHypothesis(sr7b_vdem2_sd, "autonomy1_prevv0 = autonomy1_prevv0 + d5_tt_vdem2:autonomy1_prevv0", test = "Chisq", vcov = cluster.vcov(sr7b_vdem2_sd, subset(df_global, sd_l1 == 0)$cowcode))
sr7b_vdem2_sd_equality2 <- linearHypothesis(sr7b_vdem2_sd, "autonomy1_prevv1 = autonomy1_prevv1 + d5_tt_vdem2:autonomy1_prevv1", test = "Chisq", vcov = cluster.vcov(sr7b_vdem2_sd, subset(df_global, sd_l1 == 0)$cowcode))
sr7b_vdem2_violsd_jointsig1 <- linearHypothesis(sr7b_vdem2_violsd, "autonomy1_prevv0 + d5_tt_vdem2:autonomy1_prevv0 = 0", test = "Chisq", vcov = cluster.vcov(sr7b_vdem2_violsd, subset(df_global, violsd_l1 == 0)$cowcode))
sr7b_vdem2_violsd_jointsig2 <- linearHypothesis(sr7b_vdem2_violsd, "autonomy1_prevv1 + d5_tt_vdem2:autonomy1_prevv1 = 0", test = "Chisq", vcov = cluster.vcov(sr7b_vdem2_violsd, subset(df_global, violsd_l1 == 0)$cowcode))
sr7b_vdem2_violsd_equality1 <- linearHypothesis(sr7b_vdem2_violsd, "autonomy1_prevv0 = autonomy1_prevv0 + d5_tt_vdem2:autonomy1_prevv0", test = "Chisq", vcov = cluster.vcov(sr7b_vdem2_violsd, subset(df_global, violsd_l1 == 0)$cowcode))
sr7b_vdem2_violsd_equality2 <- linearHypothesis(sr7b_vdem2_violsd, "autonomy1_prevv1 = autonomy1_prevv1 + d5_tt_vdem2:autonomy1_prevv1", test = "Chisq", vcov = cluster.vcov(sr7b_vdem2_violsd, subset(df_global, violsd_l1 == 0)$cowcode))
srg7b_vdem2_sd_jointsig1 <- linearHypothesis(srg7b_vdem2_sd, "autonomy1_prevv0 + tt_vdem2_shock:autonomy1_prevv0 = 0", test = "Chisq", vcov = cluster.vcov(srg7b_vdem2_sd, subset(df_global, sd_l1 == 0)$cowcode))
srg7b_vdem2_sd_jointsig2 <- linearHypothesis(srg7b_vdem2_sd, "autonomy1_prevv1 + tt_vdem2_shock:autonomy1_prevv1 = 0", test = "Chisq", vcov = cluster.vcov(srg7b_vdem2_sd, subset(df_global, sd_l1 == 0)$cowcode))
srg7b_vdem2_sd_equality1 <- linearHypothesis(srg7b_vdem2_sd, "autonomy1_prevv0 = autonomy1_prevv0 + tt_vdem2_shock:autonomy1_prevv0", test = "Chisq", vcov = cluster.vcov(srg7b_vdem2_sd, subset(df_global, sd_l1 == 0)$cowcode))
srg7b_vdem2_sd_equality2 <- linearHypothesis(srg7b_vdem2_sd, "autonomy1_prevv1 = autonomy1_prevv1 + tt_vdem2_shock:autonomy1_prevv1", test = "Chisq", vcov = cluster.vcov(srg7b_vdem2_sd, subset(df_global, sd_l1 == 0)$cowcode))
srg7b_vdem2_violsd_jointsig1 <- linearHypothesis(srg7b_vdem2_violsd, "autonomy1_prevv0 + tt_vdem2_shock:autonomy1_prevv0 = 0", test = "Chisq", vcov = cluster.vcov(srg7b_vdem2_violsd, subset(df_global, violsd_l1 == 0)$cowcode))
srg7b_vdem2_violsd_jointsig2 <- linearHypothesis(srg7b_vdem2_violsd, "autonomy1_prevv1 + tt_vdem2_shock:autonomy1_prevv1 = 0", test = "Chisq", vcov = cluster.vcov(srg7b_vdem2_violsd, subset(df_global, violsd_l1 == 0)$cowcode))
srg7b_vdem2_violsd_equality1 <- linearHypothesis(srg7b_vdem2_violsd, "autonomy1_prevv0 = autonomy1_prevv0 + tt_vdem2_shock:autonomy1_prevv0", test = "Chisq", vcov = cluster.vcov(srg7b_vdem2_violsd, subset(df_global, violsd_l1 == 0)$cowcode))
srg7b_vdem2_violsd_equality2 <- linearHypothesis(srg7b_vdem2_violsd, "autonomy1_prevv1 = autonomy1_prevv1 + tt_vdem2_shock:autonomy1_prevv1", test = "Chisq", vcov = cluster.vcov(srg7b_vdem2_violsd, subset(df_global, violsd_l1 == 0)$cowcode))

##########3.3 Autonomy and ethno-regional organizations (I: control for vote share of ethnic/regional parties)##########
###a) models
sr7c_vdem2_sd <- glm(as.formula(paste("sd_onset", paste(c("d5_tt_vdem2", "autonomy", "autonomy:d5_tt_vdem2", "party_reg_ethnic_vote_share", "autonomy:party_reg_ethnic_vote_share", "d5_tt_vdem2:party_reg_ethnic_vote_share","yearf","as.factor(region)",group_controls,country_controls, sdm_years), collapse = " + "), sep = " ~ ")), data=subset(df_global, sd_l1 == 0), family=binomial(link="logit"),na.action = na.exclude)
cse_sr7c_vdem2_sd <- data.frame(cluster.se(sr7c_vdem2_sd,as.factor(subset(df_global, sd_l1 == 0)$cowcode))[, 2])
sr7c_vdem2_violsd <- glm(as.formula(paste("violsd_onset", paste(c("d5_tt_vdem2", "autonomy", "autonomy:d5_tt_vdem2", "party_reg_ethnic_vote_share", "autonomy:party_reg_ethnic_vote_share", "d5_tt_vdem2:party_reg_ethnic_vote_share","yearf","as.factor(region)",group_controls,country_controls, violsdm_years), collapse = " + "), sep = " ~ ")), data=subset(df_global, violsd_l1 == 0), family=binomial(link="logit"),na.action = na.exclude)
cse_sr7c_vdem2_violsd <- data.frame(cluster.se(sr7c_vdem2_violsd,as.factor(subset(df_global, violsd_l1 == 0)$cowcode))[, 2])
srg7c_vdem2_sd <- glm(as.formula(paste("sd_onset", paste(c("tt_vdem2_shock*autonomy", "party_reg_ethnic_vote_share", "autonomy:party_reg_ethnic_vote_share", "tt_vdem2_shock:party_reg_ethnic_vote_share","yearf","as.factor(region)",group_controls,country_controls, sdm_years), collapse = " + "), sep = " ~ ")), data=subset(df_global, sd_l1 == 0), family=binomial(link="logit"),na.action = na.exclude)
cse_srg7c_vdem2_sd <- data.frame(cluster.se(srg7c_vdem2_sd,as.factor(subset(df_global, sd_l1 == 0)$cowcode))[, 2])
srg7c_vdem2_violsd <- glm(as.formula(paste("violsd_onset", paste(c("tt_vdem2_shock*autonomy", "party_reg_ethnic_vote_share", "autonomy:party_reg_ethnic_vote_share", "tt_vdem2_shock:party_reg_ethnic_vote_share","yearf","as.factor(region)",group_controls,country_controls, violsdm_years), collapse = " + "), sep = " ~ ")), data=subset(df_global, violsd_l1 == 0), family=binomial(link="logit"),na.action = na.exclude)
cse_srg7c_vdem2_violsd <- data.frame(cluster.se(srg7c_vdem2_violsd,as.factor(subset(df_global, violsd_l1 == 0)$cowcode))[, 2])
stargazer(type="text",sr7c_vdem2_sd, sr7c_vdem2_violsd, srg7c_vdem2_sd, srg7c_vdem2_violsd,se=c(cse_sr7c_vdem2_sd, cse_sr7c_vdem2_violsd, cse_srg7c_vdem2_sd, cse_srg7c_vdem2_violsd), column.labels = c("Model 1", "Model 2", "Model 3", "Model 4"),omit=c("years_no","yearf","factor"),  dep.var.labels.include=TRUE, dep.var.labels = c("SDM onset","Violent SDM onset","SDM onset","Violent SDM onset"), style = "ajps", notes.append = TRUE, notes.align = "l", model.numbers = F, notes = "Country-clustered errors in parentheses. Cubic term for time dependence of dependent variable and year- and region-fixed effects included but not reported.", order=c(paste0("^", vars_vdem2, "$"), ":", paste0("^", vars.order_main2, "$")), covariate.labels = c("Autonomy","Transition (0-5)","Autonomy x transition (0-5)","Transition proximity","Autonomy x transition proximity","Autonomy x regional party vote share","Transition (0-5) x regional party vote share","Transition proximity x regional party vote share","Most powerful","Included","Group size","TEK irredentism","TEK SDMs","Petroleum % in settlement area","Previous SDMs in country","Normalized polity score","Regional party vote share","Distance border (log)","GDP p.c. (log)","Population (log)","Ongoing SD"))
###b) Wald tests
sr7c_vdem2_sd_jointsig <- linearHypothesis(sr7c_vdem2_sd, "autonomy + d5_tt_vdem2:autonomy = 0", test = "Chisq", vcov = cluster.vcov(sr7c_vdem2_sd, subset(df_global, sd_l1 == 0)$cowcode))
sr7c_vdem2_sd_equality <- linearHypothesis(sr7c_vdem2_sd, "autonomy = autonomy + d5_tt_vdem2:autonomy", test = "Chisq", vcov = cluster.vcov(sr7c_vdem2_sd, subset(df_global, sd_l1 == 0)$cowcode))
sr7c_vdem2_violsd_jointsig <- linearHypothesis(sr7c_vdem2_violsd, "autonomy + d5_tt_vdem2:autonomy = 0", test = "Chisq", vcov = cluster.vcov(sr7c_vdem2_violsd, subset(df_global, violsd_l1 == 0)$cowcode))
sr7c_vdem2_violsd_equality <- linearHypothesis(sr7c_vdem2_violsd, "autonomy = autonomy + d5_tt_vdem2:autonomy", test = "Chisq", vcov = cluster.vcov(sr7c_vdem2_violsd, subset(df_global, violsd_l1 == 0)$cowcode))
srg7c_vdem2_sd_jointsig <- linearHypothesis(srg7c_vdem2_sd, "autonomy + tt_vdem2_shock:autonomy = 0", test = "Chisq", vcov = cluster.vcov(srg7c_vdem2_sd, subset(df_global, sd_l1 == 0)$cowcode))
srg7c_vdem2_sd_equality <- linearHypothesis(srg7c_vdem2_sd, "autonomy = autonomy + tt_vdem2_shock:autonomy", test = "Chisq", vcov = cluster.vcov(srg7c_vdem2_sd, subset(df_global, sd_l1 == 0)$cowcode))
srg7c_vdem2_violsd_jointsig <- linearHypothesis(srg7c_vdem2_violsd, "autonomy + tt_vdem2_shock:autonomy = 0", test = "Chisq", vcov = cluster.vcov(srg7c_vdem2_violsd, subset(df_global, violsd_l1 == 0)$cowcode))
srg7c_vdem2_violsd_equality <- linearHypothesis(srg7c_vdem2_violsd, "autonomy = autonomy + tt_vdem2_shock:autonomy", test = "Chisq", vcov = cluster.vcov(srg7c_vdem2_violsd, subset(df_global, violsd_l1 == 0)$cowcode))

##########3.4 Autonomy and ethno-regional organizations (II: control for number of ethnic organizations)##########
###a) models
sr7d_vdem2_sd <- glm(as.formula(paste("sd_onset", paste(c("d5_tt_vdem2", "autonomy", "autonomy:d5_tt_vdem2", "lorg_number_any_g", "autonomy:lorg_number_any_g", "d5_tt_vdem2:lorg_number_any_g","state_control","included","yearf","as.factor(region)",group_controls_org,country_controls_org, sdm_years), collapse = " + "), sep = " ~ ")), data=subset(df_global, in_epr_org == 1 & sd_l1 == 0), family=binomial(link="logit"),na.action = na.exclude)
cse_sr7d_vdem2_sd <- data.frame(cluster.se(sr7d_vdem2_sd,as.factor(subset(df_global, in_epr_org == 1 & sd_l1 == 0)$cowcode))[, 2])
sr7d_vdem2_violsd <- glm(as.formula(paste("violsd_onset", paste(c("d5_tt_vdem2", "autonomy", "autonomy:d5_tt_vdem2", "lorg_number_any_g", "autonomy:lorg_number_any_g", "d5_tt_vdem2:lorg_number_any_g","state_control","included","yearf","as.factor(region)",group_controls_org,country_controls_org, violsdm_years,"sd_l1"), collapse = " + "), sep = " ~ ")), data=subset(df_global, in_epr_org == 1 & violsd_l1 == 0), family=binomial(link="logit"),na.action = na.exclude)
cse_sr7d_vdem2_violsd <- data.frame(cluster.se(sr7d_vdem2_violsd,as.factor(subset(df_global, in_epr_org == 1 & violsd_l1 == 0)$cowcode))[, 2])
srg7d_vdem2_sd <- glm(as.formula(paste("sd_onset", paste(c("tt_vdem2_shock*autonomy", "lorg_number_any_g", "autonomy:lorg_number_any_g", "tt_vdem2_shock:lorg_number_any_g","state_control","included","yearf","as.factor(region)",group_controls_org,country_controls_org, sdm_years), collapse = " + "), sep = " ~ ")), data=subset(df_global, in_epr_org == 1 & sd_l1 == 0), family=binomial(link="logit"),na.action = na.exclude)
cse_srg7d_vdem2_sd <- data.frame(cluster.se(srg7d_vdem2_sd,as.factor(subset(df_global, in_epr_org == 1 & sd_l1 == 0)$cowcode))[, 2])
srg7d_vdem2_violsd <- glm(as.formula(paste("violsd_onset", paste(c("tt_vdem2_shock*autonomy", "lorg_number_any_g", "autonomy:lorg_number_any_g", "tt_vdem2_shock:lorg_number_any_g","state_control","included","yearf","as.factor(region)",group_controls_org,country_controls_org, violsdm_years,"sd_l1"), collapse = " + "), sep = " ~ ")), data=subset(df_global, in_epr_org == 1 & violsd_l1 == 0), family=binomial(link="logit"),na.action = na.exclude)
cse_srg7d_vdem2_violsd <- data.frame(cluster.se(srg7d_vdem2_violsd,as.factor(subset(df_global, in_epr_org == 1 & violsd_l1 == 0)$cowcode))[, 2])
stargazer(type="text",sr7d_vdem2_sd, sr7d_vdem2_violsd, srg7d_vdem2_sd, srg7d_vdem2_violsd,se=c(cse_sr7d_vdem2_sd, cse_sr7d_vdem2_violsd, cse_srg7d_vdem2_sd, cse_srg7d_vdem2_violsd), column.labels = c("Model 1", "Model 2", "Model 3", "Model 4"),omit=c("years_no","yearf","factor"),  dep.var.labels.include=TRUE, dep.var.labels = c("SDM onset","Violent SDM onset","SDM onset","Violent SDM onset"), style = "ajps", notes.append = TRUE, notes.align = "l", model.numbers = F, notes = "Country-clustered errors in parentheses. Cubic term for time dependence of dependent variable and year- and region-fixed effects included but not reported.", order=c(paste0("^", vars_vdem2, "$"), ":", paste0("^", vars.order_main2, "$")), covariate.labels = c("Autonomy","Transition (0-5)","Autonomy x transition (0-5)","Transition proximity","Autonomy x transition proximity","Autonomy x no. ethnic organizations (log)","Transition (0-5) x no. ethnic organizations (log)","Transition proximity x no. ethnic organizations (log)","Most powerful","Included","Group size","TEK irredentism","TEK SDMs","Petroleum % in settlement area","Previous SDMs in country","Normalized polity score","No. ethnic organizations (log)","Distance border (log)","GDP p.c. (log)","Population (log)","Ongoing SD"))
###b) Wald tests
sr7d_vdem2_sd_jointsig <- linearHypothesis(sr7d_vdem2_sd, "autonomy + d5_tt_vdem2:autonomy = 0", test = "Chisq", vcov = cluster.vcov(sr7d_vdem2_sd, subset(df_global, in_epr_org == 1 & sd_l1 == 0)$cowcode))
sr7d_vdem2_sd_equality <- linearHypothesis(sr7d_vdem2_sd, "autonomy = autonomy + d5_tt_vdem2:autonomy", test = "Chisq", vcov = cluster.vcov(sr7d_vdem2_sd, subset(df_global, in_epr_org == 1 & sd_l1 == 0)$cowcode))
sr7d_vdem2_violsd_jointsig <- linearHypothesis(sr7d_vdem2_violsd, "autonomy + d5_tt_vdem2:autonomy = 0", test = "Chisq", vcov = cluster.vcov(sr7d_vdem2_violsd, subset(df_global, in_epr_org == 1 & violsd_l1 == 0)$cowcode))
sr7d_vdem2_violsd_equality <- linearHypothesis(sr7d_vdem2_violsd, "autonomy = autonomy + d5_tt_vdem2:autonomy", test = "Chisq", vcov = cluster.vcov(sr7d_vdem2_violsd, subset(df_global, in_epr_org == 1 & violsd_l1 == 0)$cowcode))
srg7d_vdem2_sd_jointsig <- linearHypothesis(srg7d_vdem2_sd, "autonomy + tt_vdem2_shock:autonomy = 0", test = "Chisq", vcov = cluster.vcov(srg7d_vdem2_sd, subset(df_global, in_epr_org == 1 & sd_l1 == 0)$cowcode))
srg7d_vdem2_sd_equality <- linearHypothesis(srg7d_vdem2_sd, "autonomy = autonomy + tt_vdem2_shock:autonomy", test = "Chisq", vcov = cluster.vcov(srg7d_vdem2_sd, subset(df_global, in_epr_org == 1 & sd_l1 == 0)$cowcode))
srg7d_vdem2_violsd_jointsig <- linearHypothesis(srg7d_vdem2_violsd, "autonomy + tt_vdem2_shock:autonomy = 0", test = "Chisq", vcov = cluster.vcov(srg7d_vdem2_violsd, subset(df_global, in_epr_org == 1 & violsd_l1 == 0)$cowcode))
srg7d_vdem2_violsd_equality <- linearHypothesis(srg7d_vdem2_violsd, "autonomy = autonomy + tt_vdem2_shock:autonomy", test = "Chisq", vcov = cluster.vcov(srg7d_vdem2_violsd, subset(df_global, in_epr_org == 1 & violsd_l1 == 0)$cowcode))

##########3.5 Autonomy and economic inequality##########
###a) models
sr7e_vdem2_sd <- glm(as.formula(paste("sd_onset", paste(c("d5_tt_vdem2", "autonomy", "autonomy:d5_tt_vdem2", "low_ratio", "autonomy:low_ratio", "d5_tt_vdem2:low_ratio", "high_ratio", "autonomy:high_ratio", "d5_tt_vdem2:high_ratio","yearf","as.factor(region)",group_controls,country_controls, sdm_years), collapse = " + "), sep = " ~ ")), data=subset(df_global, sd_l1 == 0), family=binomial(link="logit"),na.action = na.exclude)
cse_sr7e_vdem2_sd <- data.frame(cluster.se(sr7e_vdem2_sd,as.factor(subset(df_global, sd_l1 == 0)$cowcode))[, 2])
sr7e_vdem2_violsd <- glm(as.formula(paste("violsd_onset", paste(c("d5_tt_vdem2", "autonomy", "autonomy:d5_tt_vdem2", "low_ratio", "autonomy:low_ratio", "d5_tt_vdem2:low_ratio", "high_ratio", "autonomy:high_ratio", "d5_tt_vdem2:high_ratio","yearf","as.factor(region)",group_controls,country_controls, violsdm_years), collapse = " + "), sep = " ~ ")), data=subset(df_global, violsd_l1 == 0), family=binomial(link="logit"),na.action = na.exclude)
cse_sr7e_vdem2_violsd <- data.frame(cluster.se(sr7e_vdem2_violsd,as.factor(subset(df_global, violsd_l1 == 0)$cowcode))[, 2])
srg7e_vdem2_sd <- glm(as.formula(paste("sd_onset", paste(c("tt_vdem2_shock*autonomy", "low_ratio", "autonomy:low_ratio", "tt_vdem2_shock:low_ratio", "high_ratio", "autonomy:high_ratio", "tt_vdem2_shock:high_ratio","yearf","as.factor(region)",group_controls,country_controls, sdm_years), collapse = " + "), sep = " ~ ")), data=subset(df_global, sd_l1 == 0), family=binomial(link="logit"),na.action = na.exclude)
cse_srg7e_vdem2_sd <- data.frame(cluster.se(srg7e_vdem2_sd,as.factor(subset(df_global, sd_l1 == 0)$cowcode))[, 2])
srg7e_vdem2_violsd <- glm(as.formula(paste("violsd_onset", paste(c("tt_vdem2_shock*autonomy", "low_ratio", "autonomy:low_ratio", "tt_vdem2_shock:low_ratio", "high_ratio", "autonomy:high_ratio", "tt_vdem2_shock:high_ratio","yearf","as.factor(region)",group_controls,country_controls, violsdm_years), collapse = " + "), sep = " ~ ")), data=subset(df_global, violsd_l1 == 0), family=binomial(link="logit"),na.action = na.exclude)
cse_srg7e_vdem2_violsd <- data.frame(cluster.se(srg7e_vdem2_violsd,as.factor(subset(df_global, violsd_l1 == 0)$cowcode))[, 2])
stargazer(type="text",sr7e_vdem2_sd, sr7e_vdem2_violsd, srg7e_vdem2_sd, srg7e_vdem2_violsd,se=c(cse_sr7e_vdem2_sd, cse_sr7e_vdem2_violsd, cse_srg7e_vdem2_sd, cse_srg7e_vdem2_violsd), column.labels = c("Model 1", "Model 2", "Model 3", "Model 4"),omit=c("years_no","yearf","factor"),  dep.var.labels.include=TRUE, dep.var.labels = c("SDM onset","Violent SDM onset","SDM onset","Violent SDM onset"), style = "ajps", notes.append = TRUE, notes.align = "l", model.numbers = F, notes = "Country-clustered errors in parentheses. Cubic term for time dependence of dependent variable and year- and region-fixed effects included but not reported.", order=c(paste0("^", vars_vdem2, "$"), ":", paste0("^", vars.order_main2, "$")), covariate.labels = c("Autonomy","Transition (0-5)","Autonomy x transition (0-5)","Transition proximity","Autonomy x transition proximity","Autonomy x low ratio","Transition (0-5) x low ratio","Transition proximity x low ratio","Autonomy x high ratio","Transition (0-5) x high ratio","Transition proximity x high ratio","Most powerful","Included","Group size","TEK irredentism","TEK SDMs","Petroleum % in settlement area","Previous SDMs in country","Normalized polity score","Low ratio","High ratio","Distance border (log)","GDP p.c. (log)","Population (log)","Ongoing SD"))
###b) Wald tests
sr7e_vdem2_sd_jointsig <- linearHypothesis(sr7e_vdem2_sd, "autonomy + d5_tt_vdem2:autonomy = 0", test = "Chisq", vcov = cluster.vcov(sr7e_vdem2_sd, subset(df_global, sd_l1 == 0)$cowcode))
sr7e_vdem2_sd_equality <- linearHypothesis(sr7e_vdem2_sd, "autonomy = autonomy + d5_tt_vdem2:autonomy", test = "Chisq", vcov = cluster.vcov(sr7e_vdem2_sd, subset(df_global, sd_l1 == 0)$cowcode))
sr7e_vdem2_violsd_jointsig <- linearHypothesis(sr7e_vdem2_violsd, "autonomy + d5_tt_vdem2:autonomy = 0", test = "Chisq", vcov = cluster.vcov(sr7e_vdem2_violsd, subset(df_global, violsd_l1 == 0)$cowcode))
sr7e_vdem2_violsd_equality <- linearHypothesis(sr7e_vdem2_violsd, "autonomy = autonomy + d5_tt_vdem2:autonomy", test = "Chisq", vcov = cluster.vcov(sr7e_vdem2_violsd, subset(df_global, violsd_l1 == 0)$cowcode))
srg7e_vdem2_sd_jointsig <- linearHypothesis(srg7e_vdem2_sd, "autonomy + tt_vdem2_shock:autonomy = 0", test = "Chisq", vcov = cluster.vcov(srg7e_vdem2_sd, subset(df_global, sd_l1 == 0)$cowcode))
srg7e_vdem2_sd_equality <- linearHypothesis(srg7e_vdem2_sd, "autonomy = autonomy + tt_vdem2_shock:autonomy", test = "Chisq", vcov = cluster.vcov(srg7e_vdem2_sd, subset(df_global, sd_l1 == 0)$cowcode))
srg7e_vdem2_violsd_jointsig <- linearHypothesis(srg7e_vdem2_violsd, "autonomy + tt_vdem2_shock:autonomy = 0", test = "Chisq", vcov = cluster.vcov(srg7e_vdem2_violsd, subset(df_global, violsd_l1 == 0)$cowcode))
srg7e_vdem2_violsd_equality <- linearHypothesis(srg7e_vdem2_violsd, "autonomy = autonomy + tt_vdem2_shock:autonomy", test = "Chisq", vcov = cluster.vcov(srg7e_vdem2_violsd, subset(df_global, violsd_l1 == 0)$cowcode))

##########3.6 Autonomy and cultural cleavages##########
###a) models
sr7f_vdem2_sd <- glm(as.formula(paste("sd_onset", paste(c("d5_tt_vdem2", "autonomy", "autonomy:d5_tt_vdem2", "cleavage_mean", "autonomy:cleavage_mean", "d5_tt_vdem2:cleavage_mean","yearf","as.factor(region)",group_controls,country_controls, sdm_years), collapse = " + "), sep = " ~ ")), data=subset(df_global, sd_l1 == 0), family=binomial(link="logit"),na.action = na.exclude)
cse_sr7f_vdem2_sd <- data.frame(cluster.se(sr7f_vdem2_sd,as.factor(subset(df_global, sd_l1 == 0)$cowcode))[, 2])
sr7f_vdem2_violsd <- glm(as.formula(paste("violsd_onset", paste(c("d5_tt_vdem2", "autonomy", "autonomy:d5_tt_vdem2", "cleavage_mean", "autonomy:cleavage_mean", "d5_tt_vdem2:cleavage_mean","yearf","as.factor(region)",group_controls,country_controls, violsdm_years), collapse = " + "), sep = " ~ ")), data=subset(df_global, violsd_l1 == 0), family=binomial(link="logit"),na.action = na.exclude)
cse_sr7f_vdem2_violsd <- data.frame(cluster.se(sr7f_vdem2_violsd,as.factor(subset(df_global, violsd_l1 == 0)$cowcode))[, 2])
srg7f_vdem2_sd <- glm(as.formula(paste("sd_onset", paste(c("tt_vdem2_shock*autonomy", "cleavage_mean", "autonomy:cleavage_mean", "tt_vdem2_shock:cleavage_mean","yearf","as.factor(region)",group_controls,country_controls, sdm_years), collapse = " + "), sep = " ~ ")), data=subset(df_global, sd_l1 == 0), family=binomial(link="logit"),na.action = na.exclude)
cse_srg7f_vdem2_sd <- data.frame(cluster.se(srg7f_vdem2_sd,as.factor(subset(df_global, sd_l1 == 0)$cowcode))[, 2])
srg7f_vdem2_violsd <- glm(as.formula(paste("violsd_onset", paste(c("tt_vdem2_shock*autonomy", "cleavage_mean", "autonomy:cleavage_mean", "tt_vdem2_shock:cleavage_mean","yearf","as.factor(region)",group_controls,country_controls, violsdm_years), collapse = " + "), sep = " ~ ")), data=subset(df_global, violsd_l1 == 0), family=binomial(link="logit"),na.action = na.exclude)
cse_srg7f_vdem2_violsd <- data.frame(cluster.se(srg7f_vdem2_violsd,as.factor(subset(df_global, violsd_l1 == 0)$cowcode))[, 2])
stargazer(type="text",sr7f_vdem2_sd, sr7f_vdem2_violsd, srg7f_vdem2_sd, srg7f_vdem2_violsd,se=c(cse_sr7f_vdem2_sd, cse_sr7f_vdem2_violsd, cse_srg7f_vdem2_sd, cse_srg7f_vdem2_violsd), column.labels = c("Model 1", "Model 2", "Model 3", "Model 4"),omit=c("years_no","yearf","factor"),  dep.var.labels.include=TRUE, dep.var.labels = c("SDM onset","Violent SDM onset","SDM onset","Violent SDM onset"), style = "ajps", notes.append = TRUE, notes.align = "l", model.numbers = F, notes = "Country-clustered errors in parentheses. Cubic term for time dependence of dependent variable and year- and region-fixed effects included but not reported.",order=c(paste0("^", vars_vdem2, "$"), ":", paste0("^", vars.order_main2, "$")),  covariate.labels = c("Autonomy","Transition (0-5)","Autonomy x transition (0-5)","Transition proximity","Autonomy x transition proximity","Autonomy x cultural cleavages","Transition (0-5) x cultural cleavages","Transition proximity x cultural cleavages","Most powerful","Included","Group size","TEK irredentism","TEK SDMs","Petroleum % in settlement area","Previous SDMs in country","Normalized polity score","Cultural cleavages","Distance border (log)","GDP p.c. (log)","Population (log)","Ongoing SD"))
###b) Wald tests
sr7f_vdem2_sd_jointsig <- linearHypothesis(sr7f_vdem2_sd, "autonomy + d5_tt_vdem2:autonomy = 0", test = "Chisq", vcov = cluster.vcov(sr7f_vdem2_sd, subset(df_global, sd_l1 == 0)$cowcode))
sr7f_vdem2_sd_equality <- linearHypothesis(sr7f_vdem2_sd, "autonomy = autonomy + d5_tt_vdem2:autonomy", test = "Chisq", vcov = cluster.vcov(sr7f_vdem2_sd, subset(df_global, sd_l1 == 0)$cowcode))
sr7f_vdem2_violsd_jointsig <- linearHypothesis(sr7f_vdem2_violsd, "autonomy + d5_tt_vdem2:autonomy = 0", test = "Chisq", vcov = cluster.vcov(sr7f_vdem2_violsd, subset(df_global, violsd_l1 == 0)$cowcode))
sr7f_vdem2_violsd_equality <- linearHypothesis(sr7f_vdem2_violsd, "autonomy = autonomy + d5_tt_vdem2:autonomy", test = "Chisq", vcov = cluster.vcov(sr7f_vdem2_violsd, subset(df_global, violsd_l1 == 0)$cowcode))
srg7f_vdem2_sd_jointsig <- linearHypothesis(srg7f_vdem2_sd, "autonomy + tt_vdem2_shock:autonomy = 0", test = "Chisq", vcov = cluster.vcov(srg7f_vdem2_sd, subset(df_global, sd_l1 == 0)$cowcode))
srg7f_vdem2_sd_equality <- linearHypothesis(srg7f_vdem2_sd, "autonomy = autonomy + tt_vdem2_shock:autonomy", test = "Chisq", vcov = cluster.vcov(srg7f_vdem2_sd, subset(df_global, sd_l1 == 0)$cowcode))
srg7f_vdem2_violsd_jointsig <- linearHypothesis(srg7f_vdem2_violsd, "autonomy + tt_vdem2_shock:autonomy = 0", test = "Chisq", vcov = cluster.vcov(srg7f_vdem2_violsd, subset(df_global, violsd_l1 == 0)$cowcode))
srg7f_vdem2_violsd_equality <- linearHypothesis(srg7f_vdem2_violsd, "autonomy = autonomy + tt_vdem2_shock:autonomy", test = "Chisq", vcov = cluster.vcov(srg7f_vdem2_violsd, subset(df_global, violsd_l1 == 0)$cowcode))

##########3.7 Autonomy and changes in government support group##########
###a) models
sr7g_vdem2_sd <- glm(as.formula(paste("sd_onset", paste(c("d5_tt_vdem2", "autonomy", "autonomy:d5_tt_vdem2", "d5_tt_govsup", "autonomy:d5_tt_govsup","included","yearf","as.factor(region)",group_controls,country_controls, sdm_years), collapse = " + "), sep = " ~ ")), data=subset(df_global, sd_l1 == 0), family=binomial(link="logit"),na.action = na.exclude)
cse_sr7g_vdem2_sd <- data.frame(cluster.se(sr7g_vdem2_sd,as.factor(subset(df_global, sd_l1 == 0)$cowcode))[, 2])
sr7g_vdem2_violsd <- glm(as.formula(paste("violsd_onset", paste(c("d5_tt_vdem2", "autonomy", "autonomy:d5_tt_vdem2", "d5_tt_govsup", "autonomy:d5_tt_govsup","included","yearf","as.factor(region)",group_controls,country_controls, violsdm_years), collapse = " + "), sep = " ~ ")), data=subset(df_global, violsd_l1 == 0), family=binomial(link="logit"),na.action = na.exclude)
cse_sr7g_vdem2_violsd <- data.frame(cluster.se(sr7g_vdem2_violsd,as.factor(subset(df_global, violsd_l1 == 0)$cowcode))[, 2])
srg7g_vdem2_sd <- glm(as.formula(paste("sd_onset", paste(c("tt_vdem2_shock*autonomy","tt_govsup_shock*autonomy","included","yearf","as.factor(region)",group_controls,country_controls, sdm_years), collapse = " + "), sep = " ~ ")), data=subset(df_global, sd_l1 == 0), family=binomial(link="logit"),na.action = na.exclude)
cse_srg7g_vdem2_sd <- data.frame(cluster.se(srg7g_vdem2_sd,as.factor(subset(df_global, sd_l1 == 0)$cowcode))[, 2])
srg7g_vdem2_violsd <- glm(as.formula(paste("violsd_onset", paste(c("tt_vdem2_shock*autonomy","tt_govsup_shock*autonomy","included","yearf","as.factor(region)",group_controls,country_controls, violsdm_years), collapse = " + "), sep = " ~ ")), data=subset(df_global, violsd_l1 == 0), family=binomial(link="logit"),na.action = na.exclude)
cse_srg7g_vdem2_violsd <- data.frame(cluster.se(srg7g_vdem2_violsd,as.factor(subset(df_global, violsd_l1 == 0)$cowcode))[, 2])
stargazer(type="text",sr7g_vdem2_sd, sr7g_vdem2_violsd, srg7g_vdem2_sd, srg7g_vdem2_violsd,se=c(cse_sr7g_vdem2_sd, cse_sr7g_vdem2_violsd, cse_srg7g_vdem2_sd, cse_srg7g_vdem2_violsd), column.labels = c("Model 1", "Model 2", "Model 3", "Model 4"),omit=c("years_no","yearf","factor"),  dep.var.labels.include=TRUE, dep.var.labels = c("SDM onset","Violent SDM onset","SDM onset","Violent SDM onset"), style = "ajps", notes.append = TRUE, notes.align = "l", model.numbers = F, notes = "Country-clustered errors in parentheses. Cubic term for time dependence of dependent variable and year- and region-fixed effects included but not reported.", order=c(paste0("^", vars_vdem2, "$"), ":", paste0("^", vars.order_main2, "$")), covariate.labels = c("Autonomy","Transition (0-5)","Autonomy x transition (0-5)","Transition proximity","Autonomy x transition proximity","Autonomy x gov. support change (0-5)","Autonomy x gov. support change proximity","Most powerful","Included","Group size","TEK irredentism","TEK SDMs","Petroleum % in settlement area","Previous SDMs in country","Normalized polity score","Gov. support change (0-5)","Gov. support change proximity","Distance border (log)","GDP p.c. (log)","Population (log)","Ongoing SD"))
###b) Wald tests
sr7g_vdem2_sd_jointsig <- linearHypothesis(sr7g_vdem2_sd, "autonomy + d5_tt_vdem2:autonomy = 0", test = "Chisq", vcov = cluster.vcov(sr7g_vdem2_sd, subset(df_global, sd_l1 == 0)$cowcode))
sr7g_vdem2_sd_equality <- linearHypothesis(sr7g_vdem2_sd, "autonomy = autonomy + d5_tt_vdem2:autonomy", test = "Chisq", vcov = cluster.vcov(sr7g_vdem2_sd, subset(df_global, sd_l1 == 0)$cowcode))
sr7g_vdem2_violsd_jointsig <- linearHypothesis(sr7g_vdem2_violsd, "autonomy + d5_tt_vdem2:autonomy = 0", test = "Chisq", vcov = cluster.vcov(sr7g_vdem2_violsd, subset(df_global, violsd_l1 == 0)$cowcode))
sr7g_vdem2_violsd_equality <- linearHypothesis(sr7g_vdem2_violsd, "autonomy = autonomy + d5_tt_vdem2:autonomy", test = "Chisq", vcov = cluster.vcov(sr7g_vdem2_violsd, subset(df_global, violsd_l1 == 0)$cowcode))
srg7g_vdem2_sd_jointsig <- linearHypothesis(srg7g_vdem2_sd, "autonomy + tt_vdem2_shock:autonomy = 0", test = "Chisq", vcov = cluster.vcov(srg7g_vdem2_sd, subset(df_global, sd_l1 == 0)$cowcode))
srg7g_vdem2_sd_equality <- linearHypothesis(srg7g_vdem2_sd, "autonomy = autonomy + tt_vdem2_shock:autonomy", test = "Chisq", vcov = cluster.vcov(srg7g_vdem2_sd, subset(df_global, sd_l1 == 0)$cowcode))
srg7g_vdem2_violsd_jointsig <- linearHypothesis(srg7g_vdem2_violsd, "autonomy + tt_vdem2_shock:autonomy = 0", test = "Chisq", vcov = cluster.vcov(srg7g_vdem2_violsd, subset(df_global, violsd_l1 == 0)$cowcode))
srg7g_vdem2_violsd_equality <- linearHypothesis(srg7g_vdem2_violsd, "autonomy = autonomy + tt_vdem2_shock:autonomy", test = "Chisq", vcov = cluster.vcov(srg7g_vdem2_violsd, subset(df_global, violsd_l1 == 0)$cowcode))


####################################################################################
#########4. Probing reverse causation and endogeneity further (appendix B.2)#########
####################################################################################

##########4.1 Excluding groups with (violent) SDM before transition##########
###a) models
sr5_vdem2_sd <- glm(as.formula(paste("sd_onset", paste(c("d5_tt_vdem2", "autonomy", "autonomy:d5_tt_vdem2","state_control","included","yearf","as.factor(region)",group_controls,country_controls, sdm_years), collapse = " + "), sep = " ~ ")), data=subset(df_global, tt_vdem2_ybt_sd2 == 0 & sd_l1 == 0), family=binomial(link="logit"),na.action = na.exclude)
cse_sr5_vdem2_sd <- data.frame(cluster.se(sr5_vdem2_sd,as.factor(subset(df_global, tt_vdem2_ybt_sd2 == 0 & sd_l1 == 0)$cowcode))[, 2])
sr5_vdem2_violsd <- glm(as.formula(paste("violsd_onset", paste(c("d5_tt_vdem2", "autonomy", "autonomy:d5_tt_vdem2","state_control","included","yearf","as.factor(region)",group_controls,country_controls, violsdm_years), collapse = " + "), sep = " ~ ")), data=subset(df_global, tt_vdem2_ybt_violsd2 == 0 & violsd_l1 == 0), family=binomial(link="logit"),na.action = na.exclude)
cse_sr5_vdem2_violsd <- data.frame(cluster.se(sr5_vdem2_violsd,as.factor(subset(df_global, tt_vdem2_ybt_violsd2 == 0 & violsd_l1 == 0)$cowcode))[, 2])
sr5b_vdem2_sd <- glm(as.formula(paste("sd_onset", paste(c("autonomy","state_control","included","yearf","as.factor(region)",group_controls,country_controls, sdm_years), collapse = " + "), sep = " ~ ")), data=subset(df_global, d5_tt_vdem2 == 1 & sd_l1 == 0 & tt_vdem2_ybt_sd2 == 0), family=binomial(link="logit"),na.action = na.exclude)
cse_sr5b_vdem2_sd <- data.frame(cluster.se(sr5b_vdem2_sd,as.factor(subset(df_global, d5_tt_vdem2 == 1 & sd_l1 == 0 & tt_vdem2_ybt_sd2 == 0)$cowcode))[, 2])
sr5b_vdem2_violsd <- glm(as.formula(paste("violsd_onset", paste(c("autonomy","state_control","included","yearf","as.factor(region)",group_controls,country_controls, violsdm_years), collapse = " + "), sep = " ~ ")), data=subset(df_global, d5_tt_vdem2 == 1 & violsd_l1 == 0 & tt_vdem2_ybt_violsd2 == 0), family=binomial(link="logit"),na.action = na.exclude)
cse_sr5b_vdem2_violsd <- data.frame(cluster.se(sr5b_vdem2_violsd,as.factor(subset(df_global, d5_tt_vdem2 == 1 & violsd_l1 == 0 & tt_vdem2_ybt_violsd2 == 0)$cowcode))[, 2])
stargazer(type="text",sr5_vdem2_sd, sr5_vdem2_violsd, sr5b_vdem2_sd, sr5b_vdem2_violsd,se=c(cse_sr5_vdem2_sd, cse_sr5_vdem2_violsd, cse_sr5b_vdem2_sd, cse_sr5b_vdem2_violsd), column.labels = c("Model 1 (all)", "Model 2 (all)", "Model 3 (transition only)", "Model 4 (transition only)"),omit=c("years_no","yearf","factor"),  dep.var.labels.include=TRUE, dep.var.labels = c("SDM onset","Violent SDM onset","SDM onset","Violent SDM onset"), style = "ajps", notes.append = TRUE, notes.align = "l", model.numbers = F, notes = "Country-clustered errors in parentheses. Cubic term for time dependence of dependent variable and year- and region-fixed effects included but not reported.", order=c(paste0("^", vars_vdem2, "$"), ":", paste0("^", vars.order_main2, "$")), covariate.labels = c("Autonomy","Transition (0-5)","Autonomy x transition (0-5)","Most powerful","Included","Group size","TEK irredentism","TEK SDMs","Petroleum % in settlement area","Previous SDMs in country","Normalized vdem2 score","Distance border (log)","GDP p.c. (log)","Population (log)","Ongoing SD"))
###b) Wald tests
sr5_vdem2_sd_jointsig <- linearHypothesis(sr5_vdem2_sd, "autonomy + d5_tt_vdem2:autonomy = 0", test = "Chisq", vcov = cluster.vcov(sr5_vdem2_sd, subset(df_global, tt_vdem2_ybt_sd2 == 0 & sd_l1 == 0)$cowcode))
sr5_vdem2_sd_equality <- linearHypothesis(sr5_vdem2_sd, "autonomy = autonomy + d5_tt_vdem2:autonomy", test = "Chisq", vcov = cluster.vcov(sr5_vdem2_sd, subset(df_global, tt_vdem2_ybt_sd2 == 0 & sd_l1 == 0)$cowcode))
sr5_vdem2_violsd_jointsig <- linearHypothesis(sr5_vdem2_violsd, "autonomy + d5_tt_vdem2:autonomy = 0", test = "Chisq", vcov = cluster.vcov(sr5_vdem2_violsd, subset(df_global, tt_vdem2_ybt_violsd2 == 0 & violsd_l1 == 0)$cowcode))
sr5_vdem2_violsd_equality <- linearHypothesis(sr5_vdem2_violsd, "autonomy = autonomy + d5_tt_vdem2:autonomy", test = "Chisq", vcov = cluster.vcov(sr5_vdem2_violsd, subset(df_global, tt_vdem2_ybt_violsd2 == 0 & violsd_l1 == 0)$cowcode))

##########4.2 Group-fixed effects##########
###a) models
lmsr3d_vdem2_sd <- lm(as.formula(paste("sd_onset", paste(c("d5_tt_vdem2", "autonomy", "autonomy:d5_tt_vdem2","state_control","included","yearf","as.factor(gwgroupid)", sdm_years), collapse = " + "), sep = " ~ ")), data=subset(df_global, sd_l1 == 0),na.action = na.exclude)
cse_lmsr3d_vdem2_sd <- data.frame(cluster.se(lmsr3d_vdem2_sd,as.factor(subset(df_global, sd_l1 == 0)$cowcode))[, 2])
lmsr3d_vdem2_violsd <- lm(as.formula(paste("violsd_onset", paste(c("d5_tt_vdem2", "autonomy", "autonomy:d5_tt_vdem2","state_control","included","yearf","as.factor(gwgroupid)", violsdm_years), collapse = " + "), sep = " ~ ")), data=subset(df_global, violsd_l1 == 0),na.action = na.exclude)
cse_lmsr3d_vdem2_violsd <- data.frame(cluster.se(lmsr3d_vdem2_violsd,as.factor(subset(df_global, violsd_l1 == 0)$cowcode))[, 2])
lmsrg3d_vdem2_sd <- lm(as.formula(paste("sd_onset", paste(c("tt_vdem2_shock*autonomy","state_control","included","yearf","as.factor(gwgroupid)", sdm_years), collapse = " + "), sep = " ~ ")), data=subset(df_global, sd_l1 == 0),na.action = na.exclude)
cse_lmsrg3d_vdem2_sd <- data.frame(cluster.se(lmsrg3d_vdem2_sd,as.factor(subset(df_global, sd_l1 == 0)$cowcode))[, 2])
lmsrg3d_vdem2_violsd <- lm(as.formula(paste("violsd_onset", paste(c("tt_vdem2_shock*autonomy","state_control","included","yearf","as.factor(gwgroupid)", violsdm_years), collapse = " + "), sep = " ~ ")), data=subset(df_global, violsd_l1 == 0),na.action = na.exclude)
cse_lmsrg3d_vdem2_violsd <- data.frame(cluster.se(lmsrg3d_vdem2_violsd,as.factor(subset(df_global, violsd_l1 == 0)$cowcode))[, 2])
stargazer(type="text",lmsr3d_vdem2_sd, lmsr3d_vdem2_violsd, lmsrg3d_vdem2_sd, lmsrg3d_vdem2_violsd,se=c(cse_lmsr3d_vdem2_sd, cse_lmsr3d_vdem2_violsd, cse_lmsrg3d_vdem2_sd, cse_lmsrg3d_vdem2_violsd), column.labels = c("Model 1", "Model 2", "Model 3", "Model 4"),omit=c("years_no","yearf","factor"),  dep.var.labels.include=TRUE, dep.var.labels = c("SDM onset","Violent SDM onset","SDM onset","Violent SDM onset"), style = "ajps", notes.append = TRUE, notes.align = "l", model.numbers = F, notes = "Country-clustered errors in parentheses. Cubic term for time dependence of dependent variable and year- and group-fixed effects included but not reported.", order=c(paste0("^", vars_vdem2, "$"), ":", paste0("^", vars.order_main2, "$")), covariate.labels = c("Autonomy","Transition (0-5)","Autonomy x transition (0-5)","Transition proximity","Autonomy x transition proximity","Most powerful","Included","Ongoing SD"))
###b) Wald tests
lmsr3d_vdem2_sd_jointsig <- linearHypothesis(lmsr3d_vdem2_sd, "autonomy + d5_tt_vdem2:autonomy = 0", test = "Chisq", vcov = cluster.vcov(lmsr3d_vdem2_sd, subset(df_global, sd_l1 == 0)$cowcode))
lmsr3d_vdem2_sd_equality <- linearHypothesis(lmsr3d_vdem2_sd, "autonomy = autonomy + d5_tt_vdem2:autonomy", test = "Chisq", vcov = cluster.vcov(lmsr3d_vdem2_sd, subset(df_global, sd_l1 == 0)$cowcode))
lmsr3d_vdem2_violsd_jointsig <- linearHypothesis(lmsr3d_vdem2_violsd, "autonomy + d5_tt_vdem2:autonomy = 0", test = "Chisq", vcov = cluster.vcov(lmsr3d_vdem2_violsd, subset(df_global, violsd_l1 == 0)$cowcode))
lmsr3d_vdem2_violsd_equality <- linearHypothesis(lmsr3d_vdem2_violsd, "autonomy = autonomy + d5_tt_vdem2:autonomy", test = "Chisq", vcov = cluster.vcov(lmsr3d_vdem2_violsd, subset(df_global, violsd_l1 == 0)$cowcode))
lmsrg3d_vdem2_sd_jointsig <- linearHypothesis(lmsrg3d_vdem2_sd, "autonomy + tt_vdem2_shock:autonomy = 0", test = "Chisq", vcov = cluster.vcov(lmsrg3d_vdem2_sd, subset(df_global, sd_l1 == 0)$cowcode))
lmsrg3d_vdem2_sd_equality <- linearHypothesis(lmsrg3d_vdem2_sd, "autonomy = autonomy + tt_vdem2_shock:autonomy", test = "Chisq", vcov = cluster.vcov(lmsrg3d_vdem2_sd, subset(df_global, sd_l1 == 0)$cowcode))
lmsrg3d_vdem2_violsd_jointsig <- linearHypothesis(lmsrg3d_vdem2_violsd, "autonomy + tt_vdem2_shock:autonomy = 0", test = "Chisq", vcov = cluster.vcov(lmsrg3d_vdem2_violsd, subset(df_global, violsd_l1 == 0)$cowcode))
lmsrg3d_vdem2_violsd_equality <- linearHypothesis(lmsrg3d_vdem2_violsd, "autonomy = autonomy + tt_vdem2_shock:autonomy", test = "Chisq", vcov = cluster.vcov(lmsrg3d_vdem2_violsd, subset(df_global, violsd_l1 == 0)$cowcode))

##########4.3 Causal sensitivity analyses##########
###a) SDM onset
sr1_vdem2_sd <- lm(as.formula(paste("sd_onset", paste(c("autonomy","included","yearf","as.factor(region)",group_controls,country_controls, sdm_years), collapse = " + "), sep = " ~ ")), data=subset(df_global, sd_l1 == 0 & d5_tt_vdem2 == 1), na.action = na.exclude)
sr1_vdem2_sd_sens <- sensemakr(sr1_vdem2_sd, "autonomy", "sd_sum_l1", kd = 1:3)
sr1_vdem2_sd_plot <- plot(sr1_vdem2_sd_sens, lim = 0.1)
###b) Violent SDM onset
sr1_vdem2_violsd <- lm(as.formula(paste("violsd_onset", paste(c("autonomy","included","yearf","as.factor(region)",group_controls,country_controls, violsdm_years), collapse = " + "), sep = " ~ ")), data=subset(df_global, violsd_l1 == 0 & d5_tt_vdem2 == 1), na.action = na.exclude)
sr1_vdem2_violsd_sens <- sensemakr(sr1_vdem2_violsd, "autonomy", "sd_sum_l1", kd = 1:3)
sr1_vdem2_violsd_plot <- plot(sr1_vdem2_violsd_sens, lim = 0.1)


####################################################################################
#########5. Unpacking transitions (appendix B.3)#########
####################################################################################

##########5.1 Differentiating between different types of transitions##########
###a) models
trans3_vdem2_sd <- glm(as.formula(paste("sd_onset", paste(c("d5_tt_vdem2_aut*autonomy","d5_tt_vdem2_dem*autonomy","d5_tt_vdem2_st*autonomy", "autonomy","state_control","included","yearf","as.factor(region)",group_controls,country_controls, sdm_years), collapse = " + "), sep = " ~ ")), data=subset(df_global, sd_l1 == 0), family=binomial(link="logit"),na.action = na.exclude)
cse_trans3_vdem2_sd <- data.frame(cluster.se(trans3_vdem2_sd,as.factor(subset(df_global, sd_l1 == 0)$cowcode))[, 2])
trans3_vdem2g_sd <- glm(as.formula(paste("sd_onset", paste(c("tt_vdem2_dem_shock*autonomy","tt_vdem2_aut_shock*autonomy","tt_vdem2_st_shock*autonomy","state_control","included","yearf","as.factor(region)",group_controls,country_controls, sdm_years), collapse = " + "), sep = " ~ ")), data=subset(df_global, sd_l1 == 0), family=binomial(link="logit"),na.action = na.exclude)
cse_trans3_vdem2g_sd <- data.frame(cluster.se(trans3_vdem2g_sd,as.factor(subset(df_global, sd_l1 == 0)$cowcode))[, 2])
trans3_vdem2_violsd <- glm(as.formula(paste("violsd_onset", paste(c("d5_tt_vdem2_aut*autonomy","d5_tt_vdem2_dem*autonomy","d5_tt_vdem2_st*autonomy", "autonomy","state_control","included","yearf","as.factor(region)",group_controls,country_controls, violsdm_years), collapse = " + "), sep = " ~ ")), data=subset(df_global, violsd_l1 == 0), family=binomial(link="logit"),na.action = na.exclude)
cse_trans3_vdem2_violsd <- data.frame(cluster.se(trans3_vdem2_violsd,as.factor(subset(df_global, violsd_l1 == 0)$cowcode))[, 2])
trans3_vdem2g_violsd <- glm(as.formula(paste("violsd_onset", paste(c("tt_vdem2_dem_shock*autonomy","tt_vdem2_aut_shock*autonomy","tt_vdem2_st_shock*autonomy","state_control","included","yearf","as.factor(region)",group_controls,country_controls, violsdm_years), collapse = " + "), sep = " ~ ")), data=subset(df_global, violsd_l1 == 0), family=binomial(link="logit"),na.action = na.exclude)
cse_trans3_vdem2g_violsd <- data.frame(cluster.se(trans3_vdem2g_violsd,as.factor(subset(df_global, violsd_l1 == 0)$cowcode))[, 2])
stargazer(type="text", trans3_vdem2_sd, trans3_vdem2_violsd, trans3_vdem2g_sd, trans3_vdem2g_violsd,se=c(cse_trans3_vdem2_sd, cse_trans3_vdem2_violsd, cse_trans3_vdem2g_sd, cse_trans3_vdem2g_violsd), column.labels = c("Model 1 (SD)","Model 1 (VIOLSD)"),omit=c("years_no","yearf","factor"),  dep.var.labels.include=TRUE, dep.var.labels = c("SDM onset","SDM onset","Violent SDM onset","Violent SDM onset"), style = "ajps", notes.append = TRUE, notes.align = "l", model.numbers = F, notes = "Country-clustered errors in parentheses. Cubic term for time dependence of dependent variable and year- and region-fixed effects included but not reported.", order=c(paste0("^", vars_vdem2_types, "$"), ":", paste0("^", vars.order_main2, "$")), covariate.labels = c("State transition (0-5)","Democratic transition (0-5)","Autocratic transition (0-5)","Autonomy","Autonomy x autocratic transition (0-5)","Autonomy x democratic transition (0-5)","Autonomy x state transition (0-5)","Autonomy x democratic transition proximity","Autonomy x autocratic transition proximity","Autonomy x state transition proximity","Most powerful","Included","Group size","TEK irredentism","TEK SDMs","Petroleum % in settlement area","Previous SDMs in country","Normalized polity score","Democratic transition proximity","Autocratic transition proximity", "State transition proximity","Distance border (log)","GDP p.c. (log)","Population (log)","Ongoing SD"))
###b) Wald tests
trans3_vdem2_sd_jointsig2 <- linearHypothesis(trans3_vdem2_sd, "autonomy + autonomy:d5_tt_vdem2_dem = 0", test = "Chisq", vcov = cluster.vcov(trans3_vdem2_sd, subset(df_global, sd_l1 == 0)$cowcode))
trans3_vdem2_sd_equality2 <- linearHypothesis(trans3_vdem2_sd, "autonomy = autonomy + autonomy:d5_tt_vdem2_dem", test = "Chisq", vcov = cluster.vcov(trans3_vdem2_sd, subset(df_global, sd_l1 == 0)$cowcode))
trans3_vdem2_sd_jointsig1 <- linearHypothesis(trans3_vdem2_sd, "autonomy + d5_tt_vdem2_aut:autonomy = 0", test = "Chisq", vcov = cluster.vcov(trans3_vdem2_sd, subset(df_global, sd_l1 == 0)$cowcode))
trans3_vdem2_sd_equality1 <- linearHypothesis(trans3_vdem2_sd, "autonomy = autonomy + d5_tt_vdem2_aut:autonomy", test = "Chisq", vcov = cluster.vcov(trans3_vdem2_sd, subset(df_global, sd_l1 == 0)$cowcode))
trans3_vdem2_sd_jointsig3 <- linearHypothesis(trans3_vdem2_sd, "autonomy + autonomy:d5_tt_vdem2_st = 0", test = "Chisq", vcov = cluster.vcov(trans3_vdem2_sd, subset(df_global, sd_l1 == 0)$cowcode))
trans3_vdem2_sd_equality3 <- linearHypothesis(trans3_vdem2_sd, "autonomy = autonomy + autonomy:d5_tt_vdem2_st", test = "Chisq", vcov = cluster.vcov(trans3_vdem2_sd, subset(df_global, sd_l1 == 0)$cowcode))
#
trans3_vdem2_violsd_jointsig2 <- linearHypothesis(trans3_vdem2_violsd, "autonomy + autonomy:d5_tt_vdem2_dem = 0", test = "Chisq", vcov = cluster.vcov(trans3_vdem2_violsd, subset(df_global, violsd_l1 == 0)$cowcode))
trans3_vdem2_violsd_equality2 <- linearHypothesis(trans3_vdem2_violsd, "autonomy = autonomy + autonomy:d5_tt_vdem2_dem", test = "Chisq", vcov = cluster.vcov(trans3_vdem2_violsd, subset(df_global, violsd_l1 == 0)$cowcode))
trans3_vdem2_violsd_jointsig1 <- linearHypothesis(trans3_vdem2_violsd, "autonomy + d5_tt_vdem2_aut:autonomy = 0", test = "Chisq", vcov = cluster.vcov(trans3_vdem2_violsd, subset(df_global, violsd_l1 == 0)$cowcode))
trans3_vdem2_violsd_equality1 <- linearHypothesis(trans3_vdem2_violsd, "autonomy = autonomy + d5_tt_vdem2_aut:autonomy", test = "Chisq", vcov = cluster.vcov(trans3_vdem2_violsd, subset(df_global, violsd_l1 == 0)$cowcode))
trans3_vdem2_violsd_jointsig3 <- linearHypothesis(trans3_vdem2_violsd, "autonomy + autonomy:d5_tt_vdem2_st = 0", test = "Chisq", vcov = cluster.vcov(trans3_vdem2_violsd, subset(df_global, violsd_l1 == 0)$cowcode))
trans3_vdem2_violsd_equality3 <- linearHypothesis(trans3_vdem2_violsd, "autonomy = autonomy + autonomy:d5_tt_vdem2_st", test = "Chisq", vcov = cluster.vcov(trans3_vdem2_violsd, subset(df_global, violsd_l1 == 0)$cowcode))
#
trans3_vdem2g_sd_jointsig2 <- linearHypothesis(trans3_vdem2g_sd, "autonomy + tt_vdem2_dem_shock:autonomy = 0", test = "Chisq", vcov = cluster.vcov(trans3_vdem2g_sd, subset(df_global, sd_l1 == 0)$cowcode))
trans3_vdem2g_sd_equality2 <- linearHypothesis(trans3_vdem2g_sd, "autonomy = autonomy + tt_vdem2_dem_shock:autonomy", test = "Chisq", vcov = cluster.vcov(trans3_vdem2g_sd, subset(df_global, sd_l1 == 0)$cowcode))
trans3_vdem2g_sd_jointsig1 <- linearHypothesis(trans3_vdem2g_sd, "autonomy + autonomy:tt_vdem2_aut_shock = 0", test = "Chisq", vcov = cluster.vcov(trans3_vdem2g_sd, subset(df_global, sd_l1 == 0)$cowcode))
trans3_vdem2g_sd_equality1 <- linearHypothesis(trans3_vdem2g_sd, "autonomy = autonomy + autonomy:tt_vdem2_aut_shock", test = "Chisq", vcov = cluster.vcov(trans3_vdem2g_sd, subset(df_global, sd_l1 == 0)$cowcode))
trans3_vdem2g_sd_jointsig3 <- linearHypothesis(trans3_vdem2g_sd, "autonomy + autonomy:tt_vdem2_st_shock = 0", test = "Chisq", vcov = cluster.vcov(trans3_vdem2g_sd, subset(df_global, sd_l1 == 0)$cowcode))
trans3_vdem2g_sd_equality3 <- linearHypothesis(trans3_vdem2g_sd, "autonomy = autonomy + autonomy:tt_vdem2_st_shock", test = "Chisq", vcov = cluster.vcov(trans3_vdem2g_sd, subset(df_global, sd_l1 == 0)$cowcode))
#
trans3_vdem2g_violsd_jointsig2 <- linearHypothesis(trans3_vdem2g_violsd, "autonomy + tt_vdem2_dem_shock:autonomy = 0", test = "Chisq", vcov = cluster.vcov(trans3_vdem2g_violsd, subset(df_global, violsd_l1 == 0)$cowcode))
trans3_vdem2g_violsd_equality2 <- linearHypothesis(trans3_vdem2g_violsd, "autonomy = autonomy + tt_vdem2_dem_shock:autonomy", test = "Chisq", vcov = cluster.vcov(trans3_vdem2g_violsd, subset(df_global, violsd_l1 == 0)$cowcode))
trans3_vdem2g_violsd_jointsig1 <- linearHypothesis(trans3_vdem2g_violsd, "autonomy + autonomy:tt_vdem2_aut_shock = 0", test = "Chisq", vcov = cluster.vcov(trans3_vdem2g_violsd, subset(df_global, violsd_l1 == 0)$cowcode))
trans3_vdem2g_violsd_equality1 <- linearHypothesis(trans3_vdem2g_violsd, "autonomy = autonomy + autonomy:tt_vdem2_aut_shock", test = "Chisq", vcov = cluster.vcov(trans3_vdem2g_violsd, subset(df_global, violsd_l1 == 0)$cowcode))
trans3_vdem2g_violsd_jointsig3 <- linearHypothesis(trans3_vdem2g_violsd, "autonomy + autonomy:tt_vdem2_st_shock = 0", test = "Chisq", vcov = cluster.vcov(trans3_vdem2g_violsd, subset(df_global, violsd_l1 == 0)$cowcode))
trans3_vdem2g_violsd_equality3 <- linearHypothesis(trans3_vdem2g_violsd, "autonomy = autonomy + autonomy:tt_vdem2_st_shock", test = "Chisq", vcov = cluster.vcov(trans3_vdem2g_violsd, subset(df_global, violsd_l1 == 0)$cowcode))

##########5.2 Alternative transition measure (based on Polity)##########
###a) models
sr1_polityc_sd <- glm(as.formula(paste("sd_onset", paste(c("d5_tt_polityc", "autonomy", "autonomy:d5_tt_polityc","state_control","included","yearf","as.factor(region)",group_controls,country_controls, sdm_years), collapse = " + "), sep = " ~ ")), data=subset(df_global, sd_l1 == 0), family=binomial(link="logit"),na.action = na.exclude)
cse_sr1_polityc_sd <- data.frame(cluster.se(sr1_polityc_sd,as.factor(subset(df_global, sd_l1 == 0)$cowcode))[, 2])
sr1_polityc_violsd <- glm(as.formula(paste("violsd_onset", paste(c("d5_tt_polityc", "autonomy", "autonomy:d5_tt_polityc","state_control","included","yearf","as.factor(region)",group_controls,country_controls, violsdm_years), collapse = " + "), sep = " ~ ")), data=subset(df_global, violsd_l1 == 0), family=binomial(link="logit"),na.action = na.exclude)
cse_sr1_polityc_violsd <- data.frame(cluster.se(sr1_polityc_violsd,as.factor(subset(df_global, violsd_l1 == 0)$cowcode))[, 2])
srg1_polityc_sd <- glm(as.formula(paste("sd_onset", paste(c("tt_polityc_shock*autonomy","state_control","included","yearf","as.factor(region)",group_controls,country_controls, sdm_years), collapse = " + "), sep = " ~ ")), data=subset(df_global, sd_l1 == 0), family=binomial(link="logit"),na.action = na.exclude)
cse_srg1_polityc_sd <- data.frame(cluster.se(srg1_polityc_sd,as.factor(subset(df_global, sd_l1 == 0)$cowcode))[, 2])
srg1_polityc_violsd <- glm(as.formula(paste("violsd_onset", paste(c("tt_polityc_shock*autonomy","state_control","included","yearf","as.factor(region)",group_controls,country_controls, violsdm_years), collapse = " + "), sep = " ~ ")), data=subset(df_global, violsd_l1 == 0), family=binomial(link="logit"),na.action = na.exclude)
cse_srg1_polityc_violsd <- data.frame(cluster.se(srg1_polityc_violsd,as.factor(subset(df_global, violsd_l1 == 0)$cowcode))[, 2])
stargazer(type="text",sr1_polityc_sd, sr1_polityc_violsd, srg1_polityc_sd, srg1_polityc_violsd,se=c(cse_sr1_polityc_sd, cse_sr1_polityc_violsd, cse_srg1_polityc_sd, cse_srg1_polityc_violsd), column.labels = c("Model 1", "Model 2", "Model 3", "Model 4"),omit=c("years_no","yearf","factor"),  dep.var.labels.include=TRUE, dep.var.labels = c("SDM onset","Violent SDM onset","SDM onset","Violent SDM onset"), style = "ajps", notes.append = TRUE, notes.align = "l", model.numbers = F, notes = "Country-clustered errors in parentheses. Cubic term for time dependence of dependent variable and year- and region-fixed effects included but not reported.", order=c(paste0("^", vars_polityc, "$"), ":", paste0("^", vars.order_main2, "$")), covariate.labels = c("Autonomy","Transition (0-5)","Autonomy x transition (0-5)","Transition proximity","Autonomy x transition proximity","Most powerful","Included","Group size","TEK irredentism","TEK SDMs","Petroleum % in settlement area","Previous SDMs in country","Normalized polity score","Distance border (log)","GDP p.c. (log)","Population (log)","Ongoing SD"))
###b) Wald tests
sr1_polityc_sd_jointsig <- linearHypothesis(sr1_polityc_sd, "autonomy + d5_tt_polityc:autonomy = 0", test = "Chisq", vcov = cluster.vcov(sr1_polityc_sd, subset(df_global, sd_l1 == 0)$cowcode))
sr1_polityc_sd_equality <- linearHypothesis(sr1_polityc_sd, "autonomy = autonomy + d5_tt_polityc:autonomy", test = "Chisq", vcov = cluster.vcov(sr1_polityc_sd, subset(df_global, sd_l1 == 0)$cowcode))
sr1_polityc_violsd_jointsig <- linearHypothesis(sr1_polityc_violsd, "autonomy + d5_tt_polityc:autonomy = 0", test = "Chisq", vcov = cluster.vcov(sr1_polityc_violsd, subset(df_global, violsd_l1 == 0)$cowcode))
sr1_polityc_violsd_equality <- linearHypothesis(sr1_polityc_violsd, "autonomy = autonomy + d5_tt_polityc:autonomy", test = "Chisq", vcov = cluster.vcov(sr1_polityc_violsd, subset(df_global, violsd_l1 == 0)$cowcode))
srg1_polityc_sd_jointsig <- linearHypothesis(srg1_polityc_sd, "autonomy + tt_polityc_shock:autonomy = 0", test = "Chisq", vcov = cluster.vcov(srg1_polityc_sd, subset(df_global, sd_l1 == 0)$cowcode))
srg1_polityc_sd_equality <- linearHypothesis(srg1_polityc_sd, "autonomy = autonomy + tt_polityc_shock:autonomy", test = "Chisq", vcov = cluster.vcov(srg1_polityc_sd, subset(df_global, sd_l1 == 0)$cowcode))
srg1_polityc_violsd_jointsig <- linearHypothesis(srg1_polityc_violsd, "autonomy + tt_polityc_shock:autonomy = 0", test = "Chisq", vcov = cluster.vcov(srg1_polityc_violsd, subset(df_global, violsd_l1 == 0)$cowcode))
srg1_polityc_violsd_equality <- linearHypothesis(srg1_polityc_violsd, "autonomy = autonomy + tt_polityc_shock:autonomy", test = "Chisq", vcov = cluster.vcov(srg1_polityc_violsd, subset(df_global, violsd_l1 == 0)$cowcode))

##########5.3 Alternative transition time windows##########
###a) models
sr2a1_vdem2_sd <- glm(as.formula(paste("sd_onset", paste(c("tt_vdem2", "autonomy", "autonomy:tt_vdem2","state_control","included","yearf","as.factor(region)",group_controls,country_controls, sdm_years), collapse = " + "), sep = " ~ ")), data=subset(df_global, sd_l1 == 0), family=binomial(link="logit"),na.action = na.exclude)
cse_sr2a1_vdem2_sd <- data.frame(cluster.se(sr2a1_vdem2_sd,as.factor(subset(df_global, sd_l1 == 0)$cowcode))[, 2])
sr2a1_vdem2_violsd <- glm(as.formula(paste("violsd_onset", paste(c("tt_vdem2", "autonomy", "autonomy:tt_vdem2","state_control","included","yearf","as.factor(region)",group_controls,country_controls, violsdm_years), collapse = " + "), sep = " ~ ")), data=subset(df_global, violsd_l1 == 0), family=binomial(link="logit"),na.action = na.exclude)
cse_sr2a1_vdem2_violsd <- data.frame(cluster.se(sr2a1_vdem2_violsd,as.factor(subset(df_global, violsd_l1 == 0)$cowcode))[, 2])
sr2a3_vdem2_sd <- glm(as.formula(paste("sd_onset", paste(c("d3_tt_vdem2", "autonomy", "autonomy:d3_tt_vdem2","state_control","included","yearf","as.factor(region)",group_controls,country_controls, sdm_years), collapse = " + "), sep = " ~ ")), data=subset(df_global, sd_l1 == 0), family=binomial(link="logit"),na.action = na.exclude)
cse_sr2a3_vdem2_sd <- data.frame(cluster.se(sr2a3_vdem2_sd,as.factor(subset(df_global, sd_l1 == 0)$cowcode))[, 2])
sr2a3_vdem2_violsd <- glm(as.formula(paste("violsd_onset", paste(c("d3_tt_vdem2", "autonomy", "autonomy:d3_tt_vdem2","state_control","included","yearf","as.factor(region)",group_controls,country_controls, violsdm_years), collapse = " + "), sep = " ~ ")), data=subset(df_global, violsd_l1 == 0), family=binomial(link="logit"),na.action = na.exclude)
cse_sr2a3_vdem2_violsd <- data.frame(cluster.se(sr2a3_vdem2_violsd,as.factor(subset(df_global, violsd_l1 == 0)$cowcode))[, 2])
stargazer(type="text", sr2a3_vdem2_sd, sr2a3_vdem2_violsd,sr2a1_vdem2_sd, sr2a1_vdem2_violsd,se=c(cse_sr2a3_vdem2_sd, cse_sr2a3_vdem2_violsd, cse_sr2a1_vdem2_sd, cse_sr2a1_vdem2_violsd), column.labels = c("Model 1", "Model 2", "Model 3", "Model 4"),omit=c("years_no","yearf","factor"),  dep.var.labels.include=TRUE, dep.var.labels = c("SDM onset","Violent SDM onset","SDM onset","Violent SDM onset"), style = "ajps", notes.append = TRUE, notes.align = "l", model.numbers = F, notes = "Country-clustered errors in parentheses. Cubic term for time dependence of dependent variable and year- and region-fixed effects included but not reported.", order=c(paste0("^", vars_vdem2_type5, "$"), ":", paste0("^", vars.order_main2, "$")), covariate.labels = c("Autonomy","Transition (0-3)","Autonomy x transition (0-3)","Transition (0)","Autonomy x transition (0)","Most powerful","Included","Group size","TEK irredentism","TEK SDMs","Petroleum % in settlement area","Previous SDMs in country","Normalized polity score","Distance border (log)","GDP p.c. (log)","Population (log)","Ongoing SD"))
##Wald tests
sr2a1_vdem2_sd_jointsig <- linearHypothesis(sr2a1_vdem2_sd, "autonomy + tt_vdem2:autonomy = 0", test = "Chisq", vcov = cluster.vcov(sr2a1_vdem2_sd, subset(df_global, sd_l1 == 0)$cowcode))
sr2a1_vdem2_sd_equality <- linearHypothesis(sr2a1_vdem2_sd, "autonomy = autonomy + tt_vdem2:autonomy", test = "Chisq", vcov = cluster.vcov(sr2a1_vdem2_sd, subset(df_global, sd_l1 == 0)$cowcode))
sr2a1_vdem2_violsd_jointsig <- linearHypothesis(sr2a1_vdem2_violsd, "autonomy + tt_vdem2:autonomy = 0", test = "Chisq", vcov = cluster.vcov(sr2a1_vdem2_violsd, subset(df_global, violsd_l1 == 0)$cowcode))
sr2a1_vdem2_violsd_equality <- linearHypothesis(sr2a1_vdem2_violsd, "autonomy = autonomy + tt_vdem2:autonomy", test = "Chisq", vcov = cluster.vcov(sr2a1_vdem2_violsd, subset(df_global, violsd_l1 == 0)$cowcode))
sr2a3_vdem2_sd_jointsig <- linearHypothesis(sr2a3_vdem2_sd, "autonomy + d3_tt_vdem2:autonomy = 0", test = "Chisq", vcov = cluster.vcov(sr2a3_vdem2_sd, subset(df_global, sd_l1 == 0)$cowcode))
sr2a3_vdem2_sd_equality <- linearHypothesis(sr2a3_vdem2_sd, "autonomy = autonomy + d3_tt_vdem2:autonomy", test = "Chisq", vcov = cluster.vcov(sr2a3_vdem2_sd, subset(df_global, sd_l1 == 0)$cowcode))
sr2a3_vdem2_violsd_jointsig <- linearHypothesis(sr2a3_vdem2_violsd, "autonomy + d3_tt_vdem2:autonomy = 0", test = "Chisq", vcov = cluster.vcov(sr2a3_vdem2_violsd, subset(df_global, violsd_l1 == 0)$cowcode))
sr2a3_vdem2_violsd_equality <- linearHypothesis(sr2a3_vdem2_violsd, "autonomy = autonomy + d3_tt_vdem2:autonomy", test = "Chisq", vcov = cluster.vcov(sr2a3_vdem2_violsd, subset(df_global, violsd_l1 == 0)$cowcode))

##########5.4 Alternative half-year values for transition proximity variable##########
###a) models
srg2b1_vdem2_sd <- glm(as.formula(paste("sd_onset", paste(c("tt_vdem2_shock1", "autonomy", "autonomy:tt_vdem2_shock1","state_control","included","yearf","as.factor(region)",group_controls,country_controls, sdm_years), collapse = " + "), sep = " ~ ")), data=subset(df_global, sd_l1 == 0), family=binomial(link="logit"),na.action = na.exclude)
cse_srg2b1_vdem2_sd <- data.frame(cluster.se(srg2b1_vdem2_sd,as.factor(subset(df_global, sd_l1 == 0)$cowcode))[, 2])
srg2b1_vdem2_violsd <- glm(as.formula(paste("violsd_onset", paste(c("tt_vdem2_shock1", "autonomy", "autonomy:tt_vdem2_shock1","state_control","included","yearf","as.factor(region)",group_controls,country_controls, violsdm_years), collapse = " + "), sep = " ~ ")), data=subset(df_global, violsd_l1 == 0), family=binomial(link="logit"),na.action = na.exclude)
cse_srg2b1_vdem2_violsd <- data.frame(cluster.se(srg2b1_vdem2_violsd,as.factor(subset(df_global, violsd_l1 == 0)$cowcode))[, 2])
srg2b3_vdem2_sd <- glm(as.formula(paste("sd_onset", paste(c("tt_vdem2_shock5", "autonomy", "autonomy:tt_vdem2_shock5","state_control","included","yearf","as.factor(region)",group_controls,country_controls, sdm_years), collapse = " + "), sep = " ~ ")), data=subset(df_global, sd_l1 == 0), family=binomial(link="logit"),na.action = na.exclude)
cse_srg2b3_vdem2_sd <- data.frame(cluster.se(srg2b3_vdem2_sd,as.factor(subset(df_global, sd_l1 == 0)$cowcode))[, 2])
srg2b3_vdem2_violsd <- glm(as.formula(paste("violsd_onset", paste(c("tt_vdem2_shock5", "autonomy", "autonomy:tt_vdem2_shock5","state_control","included","yearf","as.factor(region)",group_controls,country_controls, violsdm_years), collapse = " + "), sep = " ~ ")), data=subset(df_global, violsd_l1 == 0), family=binomial(link="logit"),na.action = na.exclude)
cse_srg2b3_vdem2_violsd <- data.frame(cluster.se(srg2b3_vdem2_violsd,as.factor(subset(df_global, violsd_l1 == 0)$cowcode))[, 2])
stargazer(type="text", srg2b3_vdem2_sd, srg2b3_vdem2_violsd,srg2b1_vdem2_sd, srg2b1_vdem2_violsd,se=c(cse_srg2b3_vdem2_sd, cse_srg2b3_vdem2_violsd, cse_srg2b1_vdem2_sd, cse_srg2b1_vdem2_violsd), column.labels = c("Model 1", "Model 2", "Model 3", "Model 4"),omit=c("years_no","yearf","factor"),  dep.var.labels.include=TRUE, dep.var.labels = c("SDM onset","Violent SDM onset","SDM onset","Violent SDM onset"), style = "ajps", notes.append = TRUE, notes.align = "l", model.numbers = F, notes = "Country-clustered errors in parentheses. Cubic term for time dependence of dependent variable and year- and region-fixed effects included but not reported.", order=c(paste0("^", vars_vdem2_shock5, "$"), ":", paste0("^", vars_vdem2_shock1, "$"), ":", paste0("^", vars.order_main2, "$")), covariate.labels = c("Autonomy","Transition proximity (d5)","Autonomy x transition proximity (d5)","Transition proximity (d1)","Autonomy x transition proximity (d1)","Most powerful","Included","Group size","TEK irredentism","TEK SDMs","Petroleum % in settlement area","Previous SDMs in country","Normalized polity score","Distance border (log)","GDP p.c. (log)","Population (log)","Ongoing SD"))
###b) Wald tests
srg2b1_vdem2_sd_jointsig <- linearHypothesis(srg2b1_vdem2_sd, "autonomy + tt_vdem2_shock1:autonomy = 0", test = "Chisq", vcov = cluster.vcov(srg2b1_vdem2_sd, subset(df_global, sd_l1 == 0)$cowcode))
srg2b1_vdem2_sd_equality <- linearHypothesis(srg2b1_vdem2_sd, "autonomy = autonomy + tt_vdem2_shock1:autonomy", test = "Chisq", vcov = cluster.vcov(srg2b1_vdem2_sd, subset(df_global, sd_l1 == 0)$cowcode))
srg2b1_vdem2_violsd_jointsig <- linearHypothesis(srg2b1_vdem2_violsd, "autonomy + tt_vdem2_shock1:autonomy = 0", test = "Chisq", vcov = cluster.vcov(srg2b1_vdem2_violsd, subset(df_global, violsd_l1 == 0)$cowcode))
srg2b1_vdem2_violsd_equality <- linearHypothesis(srg2b1_vdem2_violsd, "autonomy = autonomy + tt_vdem2_shock1:autonomy", test = "Chisq", vcov = cluster.vcov(srg2b1_vdem2_violsd, subset(df_global, violsd_l1 == 0)$cowcode))
srg2b3_vdem2_sd_jointsig <- linearHypothesis(srg2b3_vdem2_sd, "autonomy + tt_vdem2_shock5:autonomy = 0", test = "Chisq", vcov = cluster.vcov(srg2b3_vdem2_sd, subset(df_global, sd_l1 == 0)$cowcode))
srg2b3_vdem2_sd_equality <- linearHypothesis(srg2b3_vdem2_sd, "autonomy = autonomy + tt_vdem2_shock5:autonomy", test = "Chisq", vcov = cluster.vcov(srg2b3_vdem2_sd, subset(df_global, sd_l1 == 0)$cowcode))
srg2b3_vdem2_violsd_jointsig <- linearHypothesis(srg2b3_vdem2_violsd, "autonomy + tt_vdem2_shock5:autonomy = 0", test = "Chisq", vcov = cluster.vcov(srg2b3_vdem2_violsd, subset(df_global, violsd_l1 == 0)$cowcode))
srg2b3_vdem2_violsd_equality <- linearHypothesis(srg2b3_vdem2_violsd, "autonomy = autonomy + tt_vdem2_shock5:autonomy", test = "Chisq", vcov = cluster.vcov(srg2b3_vdem2_violsd, subset(df_global, violsd_l1 == 0)$cowcode))

##########5.5 Control for anocratic regime type##########
###a) models
sr4a_vdem2_sd <- glm(as.formula(paste("sd_onset", paste(c("d5_tt_vdem2", "autonomy", "autonomy:d5_tt_vdem2","state_control","included","yearf","as.factor(region)","I(polityca^2)",group_controls,country_controls, sdm_years), collapse = " + "), sep = " ~ ")), data=subset(df_global, sd_l1 == 0), family=binomial(link="logit"),na.action = na.exclude)
cse_sr4a_vdem2_sd <- data.frame(cluster.se(sr4a_vdem2_sd,as.factor(subset(df_global, sd_l1 == 0)$cowcode))[, 2])
sr4a_vdem2_violsd <- glm(as.formula(paste("violsd_onset", paste(c("d5_tt_vdem2", "autonomy", "autonomy:d5_tt_vdem2","state_control","included","yearf","as.factor(region)","I(polityca^2)",group_controls,country_controls, violsdm_years), collapse = " + "), sep = " ~ ")), data=subset(df_global, violsd_l1 == 0), family=binomial(link="logit"),na.action = na.exclude)
cse_sr4a_vdem2_violsd <- data.frame(cluster.se(sr4a_vdem2_violsd,as.factor(subset(df_global, violsd_l1 == 0)$cowcode))[, 2])
srg4a_vdem2_sd <- glm(as.formula(paste("sd_onset", paste(c("tt_vdem2_shock*autonomy","state_control","included","yearf","as.factor(region)","I(polityca^2)",group_controls,country_controls, sdm_years), collapse = " + "), sep = " ~ ")), data=subset(df_global, sd_l1 == 0), family=binomial(link="logit"),na.action = na.exclude)
cse_srg4a_vdem2_sd <- data.frame(cluster.se(srg4a_vdem2_sd,as.factor(subset(df_global, sd_l1 == 0)$cowcode))[, 2])
srg4a_vdem2_violsd <- glm(as.formula(paste("violsd_onset", paste(c("tt_vdem2_shock*autonomy","state_control","included","yearf","as.factor(region)","I(polityca^2)",group_controls,country_controls, violsdm_years), collapse = " + "), sep = " ~ ")), data=subset(df_global, violsd_l1 == 0), family=binomial(link="logit"),na.action = na.exclude)
cse_srg4a_vdem2_violsd <- data.frame(cluster.se(srg4a_vdem2_violsd,as.factor(subset(df_global, violsd_l1 == 0)$cowcode))[, 2])
stargazer(type="text",sr4a_vdem2_sd, sr4a_vdem2_violsd, srg4a_vdem2_sd, srg4a_vdem2_violsd,se=c(cse_sr4a_vdem2_sd, cse_sr4a_vdem2_violsd, cse_srg4a_vdem2_sd, cse_srg4a_vdem2_violsd), column.labels = c("Model 1", "Model 2", "Model 3", "Model 4"),omit=c("years_no","yearf","factor"),  dep.var.labels.include=TRUE, dep.var.labels = c("SDM onset","Violent SDM onset","SDM onset","Violent SDM onset"), style = "ajps", notes.append = TRUE, notes.align = "l", model.numbers = F, notes = "Country-clustered errors in parentheses. Cubic term for time dependence of dependent variable and year- and region-fixed effects included but not reported.", order=c(paste0("^", vars_vdem2, "$"), ":", paste0("^", vars.order_main2, "$")), covariate.labels = c("Autonomy","Transition (0-5)","Autonomy x transition (0-5)","Transition proximity","Autonomy x transition proximity","Most powerful","Included","Group size","TEK irredentism","TEK SDMs","Petroleum % in settlement area","Previous SDMs in country","Normalized polity score","Normalized polity score 2","Distance border (log)","GDP p.c. (log)","Population (log)","Ongoing SD"))
##Wald tests
sr4a_vdem2_sd_jointsig <- linearHypothesis(sr4a_vdem2_sd, "autonomy + d5_tt_vdem2:autonomy = 0", test = "Chisq", vcov = cluster.vcov(sr4a_vdem2_sd, subset(df_global, sd_l1 == 0)$cowcode))
sr4a_vdem2_sd_equality <- linearHypothesis(sr4a_vdem2_sd, "autonomy = autonomy + d5_tt_vdem2:autonomy", test = "Chisq", vcov = cluster.vcov(sr4a_vdem2_sd, subset(df_global, sd_l1 == 0)$cowcode))
sr4a_vdem2_violsd_jointsig <- linearHypothesis(sr4a_vdem2_violsd, "autonomy + d5_tt_vdem2:autonomy = 0", test = "Chisq", vcov = cluster.vcov(sr4a_vdem2_violsd, subset(df_global, violsd_l1 == 0)$cowcode))
sr4a_vdem2_violsd_equality <- linearHypothesis(sr4a_vdem2_violsd, "autonomy = autonomy + d5_tt_vdem2:autonomy", test = "Chisq", vcov = cluster.vcov(sr4a_vdem2_violsd, subset(df_global, violsd_l1 == 0)$cowcode))
srg4a_vdem2_sd_jointsig <- linearHypothesis(srg4a_vdem2_sd, "autonomy + tt_vdem2_shock:autonomy = 0", test = "Chisq", vcov = cluster.vcov(srg4a_vdem2_sd, subset(df_global, sd_l1 == 0)$cowcode))
srg4a_vdem2_sd_equality <- linearHypothesis(srg4a_vdem2_sd, "autonomy = autonomy + tt_vdem2_shock:autonomy", test = "Chisq", vcov = cluster.vcov(srg4a_vdem2_sd, subset(df_global, sd_l1 == 0)$cowcode))
srg4a_vdem2_violsd_jointsig <- linearHypothesis(srg4a_vdem2_violsd, "autonomy + tt_vdem2_shock:autonomy = 0", test = "Chisq", vcov = cluster.vcov(srg4a_vdem2_violsd, subset(df_global, violsd_l1 == 0)$cowcode))
srg4a_vdem2_violsd_equality <- linearHypothesis(srg4a_vdem2_violsd, "autonomy = autonomy + tt_vdem2_shock:autonomy", test = "Chisq", vcov = cluster.vcov(srg4a_vdem2_violsd, subset(df_global, violsd_l1 == 0)$cowcode))


####################################################################################
#########6. Alteratios to the dependent variable, sample, and specification (appendix B.4)#########
####################################################################################

##########6.1 SDM radicalization##########
###a) models
sr1_vdem2_inc_sd <- glm(as.formula(paste("inc_sd_onset", paste(c("d5_tt_vdem2", "autonomy", "autonomy:d5_tt_vdem2","state_control","included","yearf","as.factor(region)",group_controls,country_controls, sdm_years,"sd_l1"), collapse = " + "), sep = " ~ ")), data=subset(df_global, irre_sd_l1 == 0), family=binomial(link="logit"),na.action = na.exclude)
cse_sr1_vdem2_inc_sd <- data.frame(cluster.se(sr1_vdem2_inc_sd,as.factor(subset(df_global, irre_sd_l1 == 0)$cowcode))[, 2])
srg1_vdem2_inc_sd <- glm(as.formula(paste("inc_sd_onset", paste(c("tt_vdem2_shock*autonomy","state_control","included","yearf","as.factor(region)",group_controls,country_controls, sdm_years,"sd_l1"), collapse = " + "), sep = " ~ ")), data=subset(df_global, irre_sd_l1 == 0), family=binomial(link="logit"),na.action = na.exclude)
cse_srg1_vdem2_inc_sd <- data.frame(cluster.se(srg1_vdem2_inc_sd,as.factor(subset(df_global, irre_sd_l1 == 0)$cowcode))[, 2])
stargazer(type="text",sr1_vdem2_inc_sd, srg1_vdem2_inc_sd,se=c(cse_sr1_vdem2_inc_sd, cse_srg1_vdem2_inc_sd), column.labels = c("Model 1", "Model 2", "Model 3", "Model 4"),omit=c("years_no","yearf","factor"),  dep.var.labels.include=TRUE, dep.var.labels = c("SDM onset","Violent SDM onset","SDM onset","Violent SDM onset"), style = "ajps", notes.append = TRUE, notes.align = "l", model.numbers = F, notes = "Country-clustered errors in parentheses. Cubic term for time dependence of dependent variable and year- and region-fixed effects included but not reported.", order=c(paste0("^", vars_vdem2, "$"), ":", paste0("^", vars.order_main2, "$")), covariate.labels = c("Autonomy","Transition (0-5)","Autonomy x transition (0-5)","Transition proximity","Autonomy x transition proximity","Most powerful","Included","Group size","TEK irredentism","TEK SDMs","Petroleum % in settlement area","Previous SDMs in country","Normalized polity score","Distance border (log)","GDP p.c. (log)","Population (log)","Ongoing SD"))
##Wald tests
sr1_vdem2_inc_sd_jointsig <- linearHypothesis(sr1_vdem2_inc_sd, "autonomy + d5_tt_vdem2:autonomy = 0", test = "Chisq", vcov = cluster.vcov(sr1_vdem2_inc_sd, subset(df_global, irre_sd_l1 == 0)$cowcode))
sr1_vdem2_inc_sd_equality <- linearHypothesis(sr1_vdem2_inc_sd, "autonomy = autonomy + d5_tt_vdem2:autonomy", test = "Chisq", vcov = cluster.vcov(sr1_vdem2_inc_sd, subset(df_global, irre_sd_l1 == 0)$cowcode))
srg1_vdem2_inc_sd_jointsig <- linearHypothesis(srg1_vdem2_inc_sd, "autonomy + tt_vdem2_shock:autonomy = 0", test = "Chisq", vcov = cluster.vcov(srg1_vdem2_inc_sd, subset(df_global, irre_sd_l1 == 0)$cowcode))
srg1_vdem2_inc_sd_equality <- linearHypothesis(srg1_vdem2_inc_sd, "autonomy = autonomy + tt_vdem2_shock:autonomy", test = "Chisq", vcov = cluster.vcov(srg1_vdem2_inc_sd, subset(df_global, irre_sd_l1 == 0)$cowcode))

##########6.2 Excluding CEE-FSU post-1989 transition##########
###a) models
sr3b_vdem2_sd <- glm(as.formula(paste("sd_onset", paste(c("d5_tt_vdem2", "autonomy", "autonomy:d5_tt_vdem2","state_control","included","yearf","as.factor(region)",group_controls,country_controls, sdm_years), collapse = " + "), sep = " ~ ")), data=subset(df_global, sd_l1 == 0 & !((cowcode == 290 | cowcode == 310 | cowcode == 315 | cowcode == 316 | cowcode == 317 | cowcode == 339 | cowcode == 341 | cowcode == 343 | cowcode == 344 | cowcode == 345 | cowcode == 346 | cowcode == 347 | cowcode == 349 | cowcode == 355 | cowcode == 359 | cowcode == 360 | cowcode == 365 | cowcode == 366 | cowcode == 367 | cowcode == 368 | cowcode == 369 | cowcode == 370 | cowcode == 371 | cowcode == 372 | cowcode == 373 | cowcode == 701 | cowcode == 702 | cowcode == 703 | cowcode == 704 | cowcode == 705)  & year >= 1990)), family=binomial(link="logit"),na.action = na.exclude)
cse_sr3b_vdem2_sd <- data.frame(cluster.se(sr3b_vdem2_sd,as.factor(subset(df_global, sd_l1 == 0 & !((cowcode == 290 | cowcode == 310 | cowcode == 315 | cowcode == 316 | cowcode == 317 | cowcode == 339 | cowcode == 341 | cowcode == 343 | cowcode == 344 | cowcode == 345 | cowcode == 346 | cowcode == 347 | cowcode == 349 | cowcode == 355 | cowcode == 359 | cowcode == 360 | cowcode == 365 | cowcode == 366 | cowcode == 367 | cowcode == 368 | cowcode == 369 | cowcode == 370 | cowcode == 371 | cowcode == 372 | cowcode == 373 | cowcode == 701 | cowcode == 702 | cowcode == 703 | cowcode == 704 | cowcode == 705)  & year >= 1990))$cowcode))[, 2])
sr3b_vdem2_violsd <- glm(as.formula(paste("violsd_onset", paste(c("d5_tt_vdem2", "autonomy", "autonomy:d5_tt_vdem2","state_control","included","yearf","as.factor(region)",group_controls,country_controls, violsdm_years), collapse = " + "), sep = " ~ ")), data=subset(df_global, violsd_l1 == 0 & !((cowcode == 290 | cowcode == 310 | cowcode == 315 | cowcode == 316 | cowcode == 317 | cowcode == 339 | cowcode == 341 | cowcode == 343 | cowcode == 344 | cowcode == 345 | cowcode == 346 | cowcode == 347 | cowcode == 349 | cowcode == 355 | cowcode == 359 | cowcode == 360 | cowcode == 365 | cowcode == 366 | cowcode == 367 | cowcode == 368 | cowcode == 369 | cowcode == 370 | cowcode == 371 | cowcode == 372 | cowcode == 373 | cowcode == 701 | cowcode == 702 | cowcode == 703 | cowcode == 704 | cowcode == 705)  & year >= 1990)), family=binomial(link="logit"),na.action = na.exclude)
cse_sr3b_vdem2_violsd <- data.frame(cluster.se(sr3b_vdem2_violsd,as.factor(subset(df_global, violsd_l1 == 0 & !((cowcode == 290 | cowcode == 310 | cowcode == 315 | cowcode == 316 | cowcode == 317 | cowcode == 339 | cowcode == 341 | cowcode == 343 | cowcode == 344 | cowcode == 345 | cowcode == 346 | cowcode == 347 | cowcode == 349 | cowcode == 355 | cowcode == 359 | cowcode == 360 | cowcode == 365 | cowcode == 366 | cowcode == 367 | cowcode == 368 | cowcode == 369 | cowcode == 370 | cowcode == 371 | cowcode == 372 | cowcode == 373 | cowcode == 701 | cowcode == 702 | cowcode == 703 | cowcode == 704 | cowcode == 705)  & year >= 1990))$cowcode))[, 2])
srg3b_vdem2_sd <- glm(as.formula(paste("sd_onset", paste(c("tt_vdem2_shock*autonomy","state_control","included","yearf","as.factor(region)",group_controls,country_controls, sdm_years), collapse = " + "), sep = " ~ ")), data=subset(df_global, sd_l1 == 0 & !((cowcode == 290 | cowcode == 310 | cowcode == 315 | cowcode == 316 | cowcode == 317 | cowcode == 339 | cowcode == 341 | cowcode == 343 | cowcode == 344 | cowcode == 345 | cowcode == 346 | cowcode == 347 | cowcode == 349 | cowcode == 355 | cowcode == 359 | cowcode == 360 | cowcode == 365 | cowcode == 366 | cowcode == 367 | cowcode == 368 | cowcode == 369 | cowcode == 370 | cowcode == 371 | cowcode == 372 | cowcode == 373 | cowcode == 701 | cowcode == 702 | cowcode == 703 | cowcode == 704 | cowcode == 705)  & year >= 1990)), family=binomial(link="logit"),na.action = na.exclude)
cse_srg3b_vdem2_sd <- data.frame(cluster.se(srg3b_vdem2_sd,as.factor(subset(df_global, sd_l1 == 0 & !((cowcode == 290 | cowcode == 310 | cowcode == 315 | cowcode == 316 | cowcode == 317 | cowcode == 339 | cowcode == 341 | cowcode == 343 | cowcode == 344 | cowcode == 345 | cowcode == 346 | cowcode == 347 | cowcode == 349 | cowcode == 355 | cowcode == 359 | cowcode == 360 | cowcode == 365 | cowcode == 366 | cowcode == 367 | cowcode == 368 | cowcode == 369 | cowcode == 370 | cowcode == 371 | cowcode == 372 | cowcode == 373 | cowcode == 701 | cowcode == 702 | cowcode == 703 | cowcode == 704 | cowcode == 705)  & year >= 1990))$cowcode))[, 2])
srg3b_vdem2_violsd <- glm(as.formula(paste("violsd_onset", paste(c("tt_vdem2_shock*autonomy","state_control","included","yearf","as.factor(region)",group_controls,country_controls, violsdm_years), collapse = " + "), sep = " ~ ")), data=subset(df_global, violsd_l1 == 0 & !((cowcode == 290 | cowcode == 310 | cowcode == 315 | cowcode == 316 | cowcode == 317 | cowcode == 339 | cowcode == 341 | cowcode == 343 | cowcode == 344 | cowcode == 345 | cowcode == 346 | cowcode == 347 | cowcode == 349 | cowcode == 355 | cowcode == 359 | cowcode == 360 | cowcode == 365 | cowcode == 366 | cowcode == 367 | cowcode == 368 | cowcode == 369 | cowcode == 370 | cowcode == 371 | cowcode == 372 | cowcode == 373 | cowcode == 701 | cowcode == 702 | cowcode == 703 | cowcode == 704 | cowcode == 705)  & year >= 1990)), family=binomial(link="logit"),na.action = na.exclude)
cse_srg3b_vdem2_violsd <- data.frame(cluster.se(srg3b_vdem2_violsd,as.factor(subset(df_global, violsd_l1 == 0 & !((cowcode == 290 | cowcode == 310 | cowcode == 315 | cowcode == 316 | cowcode == 317 | cowcode == 339 | cowcode == 341 | cowcode == 343 | cowcode == 344 | cowcode == 345 | cowcode == 346 | cowcode == 347 | cowcode == 349 | cowcode == 355 | cowcode == 359 | cowcode == 360 | cowcode == 365 | cowcode == 366 | cowcode == 367 | cowcode == 368 | cowcode == 369 | cowcode == 370 | cowcode == 371 | cowcode == 372 | cowcode == 373 | cowcode == 701 | cowcode == 702 | cowcode == 703 | cowcode == 704 | cowcode == 705)  & year >= 1990))$cowcode))[, 2])
stargazer(type="text",sr3b_vdem2_sd, sr3b_vdem2_violsd, srg3b_vdem2_sd, srg3b_vdem2_violsd,se=c(cse_sr3b_vdem2_sd, cse_sr3b_vdem2_violsd, cse_srg3b_vdem2_sd, cse_srg3b_vdem2_violsd), column.labels = c("Model 1", "Model 2", "Model 3", "Model 4"),omit=c("years_no","yearf","factor"),  dep.var.labels.include=TRUE, dep.var.labels = c("SDM onset","Violent SDM onset","SDM onset","Violent SDM onset"), style = "ajps", notes.append = TRUE, notes.align = "l", model.numbers = F, notes = "Country-clustered errors in parentheses. Cubic term for time dependence of dependent variable and year- and region-fixed effects included but not reported.", order=c(paste0("^", vars_vdem2, "$"), ":", paste0("^", vars.order_main2, "$")), covariate.labels = c("Autonomy","Transition (0-5)","Autonomy x transition (0-5)","Transition proximity","Autonomy x transition proximity","Most powerful","Included","Group size","TEK irredentism","TEK SDMs","Petroleum % in settlement area","Previous SDMs in country","Normalized polity score","Distance border (log)","GDP p.c. (log)","Population (log)","Ongoing SD"))
###b) Wald tests
sr3b_vdem2_sd_jointsig <- linearHypothesis(sr3b_vdem2_sd, "autonomy + d5_tt_vdem2:autonomy = 0", test = "Chisq", vcov = cluster.vcov(sr3b_vdem2_sd, subset(df_global, sd_l1 == 0& !((cowcode == 290 | cowcode == 310 | cowcode == 315 | cowcode == 316 | cowcode == 317 | cowcode == 339 | cowcode == 341 | cowcode == 343 | cowcode == 344 | cowcode == 345 | cowcode == 346 | cowcode == 347 | cowcode == 349 | cowcode == 355 | cowcode == 359 | cowcode == 360 | cowcode == 365 | cowcode == 366 | cowcode == 367 | cowcode == 368 | cowcode == 369 | cowcode == 370 | cowcode == 371 | cowcode == 372 | cowcode == 373 | cowcode == 701 | cowcode == 702 | cowcode == 703 | cowcode == 704 | cowcode == 705)  & year >= 1990))$cowcode))
sr3b_vdem2_sd_equality <- linearHypothesis(sr3b_vdem2_sd, "autonomy = autonomy + d5_tt_vdem2:autonomy", test = "Chisq", vcov = cluster.vcov(sr3b_vdem2_sd, subset(df_global, sd_l1 == 0& !((cowcode == 290 | cowcode == 310 | cowcode == 315 | cowcode == 316 | cowcode == 317 | cowcode == 339 | cowcode == 341 | cowcode == 343 | cowcode == 344 | cowcode == 345 | cowcode == 346 | cowcode == 347 | cowcode == 349 | cowcode == 355 | cowcode == 359 | cowcode == 360 | cowcode == 365 | cowcode == 366 | cowcode == 367 | cowcode == 368 | cowcode == 369 | cowcode == 370 | cowcode == 371 | cowcode == 372 | cowcode == 373 | cowcode == 701 | cowcode == 702 | cowcode == 703 | cowcode == 704 | cowcode == 705)  & year >= 1990))$cowcode))
sr3b_vdem2_violsd_jointsig <- linearHypothesis(sr3b_vdem2_violsd, "autonomy + d5_tt_vdem2:autonomy = 0", test = "Chisq", vcov = cluster.vcov(sr3b_vdem2_violsd, subset(df_global, violsd_l1 == 0& !((cowcode == 290 | cowcode == 310 | cowcode == 315 | cowcode == 316 | cowcode == 317 | cowcode == 339 | cowcode == 341 | cowcode == 343 | cowcode == 344 | cowcode == 345 | cowcode == 346 | cowcode == 347 | cowcode == 349 | cowcode == 355 | cowcode == 359 | cowcode == 360 | cowcode == 365 | cowcode == 366 | cowcode == 367 | cowcode == 368 | cowcode == 369 | cowcode == 370 | cowcode == 371 | cowcode == 372 | cowcode == 373 | cowcode == 701 | cowcode == 702 | cowcode == 703 | cowcode == 704 | cowcode == 705)  & year >= 1990))$cowcode))
sr3b_vdem2_violsd_equality <- linearHypothesis(sr3b_vdem2_violsd, "autonomy = autonomy + d5_tt_vdem2:autonomy", test = "Chisq", vcov = cluster.vcov(sr3b_vdem2_violsd, subset(df_global, violsd_l1 == 0& !((cowcode == 290 | cowcode == 310 | cowcode == 315 | cowcode == 316 | cowcode == 317 | cowcode == 339 | cowcode == 341 | cowcode == 343 | cowcode == 344 | cowcode == 345 | cowcode == 346 | cowcode == 347 | cowcode == 349 | cowcode == 355 | cowcode == 359 | cowcode == 360 | cowcode == 365 | cowcode == 366 | cowcode == 367 | cowcode == 368 | cowcode == 369 | cowcode == 370 | cowcode == 371 | cowcode == 372 | cowcode == 373 | cowcode == 701 | cowcode == 702 | cowcode == 703 | cowcode == 704 | cowcode == 705)  & year >= 1990))$cowcode))
srg3b_vdem2_sd_jointsig <- linearHypothesis(srg3b_vdem2_sd, "autonomy + tt_vdem2_shock:autonomy = 0", test = "Chisq", vcov = cluster.vcov(srg3b_vdem2_sd, subset(df_global, sd_l1 == 0& !((cowcode == 290 | cowcode == 310 | cowcode == 315 | cowcode == 316 | cowcode == 317 | cowcode == 339 | cowcode == 341 | cowcode == 343 | cowcode == 344 | cowcode == 345 | cowcode == 346 | cowcode == 347 | cowcode == 349 | cowcode == 355 | cowcode == 359 | cowcode == 360 | cowcode == 365 | cowcode == 366 | cowcode == 367 | cowcode == 368 | cowcode == 369 | cowcode == 370 | cowcode == 371 | cowcode == 372 | cowcode == 373 | cowcode == 701 | cowcode == 702 | cowcode == 703 | cowcode == 704 | cowcode == 705)  & year >= 1990))$cowcode))
srg3b_vdem2_sd_equality <- linearHypothesis(srg3b_vdem2_sd, "autonomy = autonomy + tt_vdem2_shock:autonomy", test = "Chisq", vcov = cluster.vcov(srg3b_vdem2_sd, subset(df_global, sd_l1 == 0& !((cowcode == 290 | cowcode == 310 | cowcode == 315 | cowcode == 316 | cowcode == 317 | cowcode == 339 | cowcode == 341 | cowcode == 343 | cowcode == 344 | cowcode == 345 | cowcode == 346 | cowcode == 347 | cowcode == 349 | cowcode == 355 | cowcode == 359 | cowcode == 360 | cowcode == 365 | cowcode == 366 | cowcode == 367 | cowcode == 368 | cowcode == 369 | cowcode == 370 | cowcode == 371 | cowcode == 372 | cowcode == 373 | cowcode == 701 | cowcode == 702 | cowcode == 703 | cowcode == 704 | cowcode == 705)  & year >= 1990))$cowcode))
srg3b_vdem2_violsd_jointsig <- linearHypothesis(srg3b_vdem2_violsd, "autonomy + tt_vdem2_shock:autonomy = 0", test = "Chisq", vcov = cluster.vcov(srg3b_vdem2_violsd, subset(df_global, violsd_l1 == 0& !((cowcode == 290 | cowcode == 310 | cowcode == 315 | cowcode == 316 | cowcode == 317 | cowcode == 339 | cowcode == 341 | cowcode == 343 | cowcode == 344 | cowcode == 345 | cowcode == 346 | cowcode == 347 | cowcode == 349 | cowcode == 355 | cowcode == 359 | cowcode == 360 | cowcode == 365 | cowcode == 366 | cowcode == 367 | cowcode == 368 | cowcode == 369 | cowcode == 370 | cowcode == 371 | cowcode == 372 | cowcode == 373 | cowcode == 701 | cowcode == 702 | cowcode == 703 | cowcode == 704 | cowcode == 705)  & year >= 1990))$cowcode))
srg3b_vdem2_violsd_equality <- linearHypothesis(srg3b_vdem2_violsd, "autonomy = autonomy + tt_vdem2_shock:autonomy", test = "Chisq", vcov = cluster.vcov(srg3b_vdem2_violsd, subset(df_global, violsd_l1 == 0& !((cowcode == 290 | cowcode == 310 | cowcode == 315 | cowcode == 316 | cowcode == 317 | cowcode == 339 | cowcode == 341 | cowcode == 343 | cowcode == 344 | cowcode == 345 | cowcode == 346 | cowcode == 347 | cowcode == 349 | cowcode == 355 | cowcode == 359 | cowcode == 360 | cowcode == 365 | cowcode == 366 | cowcode == 367 | cowcode == 368 | cowcode == 369 | cowcode == 370 | cowcode == 371 | cowcode == 372 | cowcode == 373 | cowcode == 701 | cowcode == 702 | cowcode == 703 | cowcode == 704 | cowcode == 705)  & year >= 1990))$cowcode))

##########6.2 Excluding new states##########
###a) models
sr3c_vdem2_st_sd <- glm(as.formula(paste("sd_onset", paste(c("d5_tt_vdem2", "autonomy", "autonomy:d5_tt_vdem2","state_control","included","yearf","as.factor(region)",group_controls,country_controls, sdm_years), collapse = " + "), sep = " ~ ")), data=subset(df_global, sd_l1 == 0 & d5_tt_vdem2_st == 0), family=binomial(link="logit"),na.action = na.exclude)
cse_sr3c_vdem2_st_sd <- data.frame(cluster.se(sr3c_vdem2_st_sd,as.factor(subset(df_global, sd_l1 == 0 & d5_tt_vdem2_st == 0)$cowcode))[, 2])
sr3c_vdem2_st_violsd <- glm(as.formula(paste("violsd_onset", paste(c("d5_tt_vdem2", "autonomy", "autonomy:d5_tt_vdem2","state_control","included","yearf","as.factor(region)",group_controls,country_controls, violsdm_years), collapse = " + "), sep = " ~ ")), data=subset(df_global, violsd_l1 == 0 & d5_tt_vdem2_st == 0), family=binomial(link="logit"),na.action = na.exclude)
cse_sr3c_vdem2_st_violsd <- data.frame(cluster.se(sr3c_vdem2_st_violsd,as.factor(subset(df_global, violsd_l1 == 0 & d5_tt_vdem2_st == 0)$cowcode))[, 2])
srg3c_vdem2_st_sd <- glm(as.formula(paste("sd_onset", paste(c("tt_vdem2_shock*autonomy","state_control","included","yearf","as.factor(region)",group_controls,country_controls, sdm_years), collapse = " + "), sep = " ~ ")), data=subset(df_global, sd_l1 == 0 & d5_tt_vdem2_st == 0), family=binomial(link="logit"),na.action = na.exclude)
cse_srg3c_vdem2_st_sd <- data.frame(cluster.se(srg3c_vdem2_st_sd,as.factor(subset(df_global, sd_l1 == 0 & d5_tt_vdem2_st == 0)$cowcode))[, 2])
srg3c_vdem2_st_violsd <- glm(as.formula(paste("violsd_onset", paste(c("tt_vdem2_shock*autonomy","state_control","included","yearf","as.factor(region)",group_controls,country_controls, violsdm_years), collapse = " + "), sep = " ~ ")), data=subset(df_global, violsd_l1 == 0 & d5_tt_vdem2_st == 0), family=binomial(link="logit"),na.action = na.exclude)
cse_srg3c_vdem2_st_violsd <- data.frame(cluster.se(srg3c_vdem2_st_violsd,as.factor(subset(df_global, violsd_l1 == 0 & d5_tt_vdem2_st == 0)$cowcode))[, 2])
stargazer(type="text",sr3c_vdem2_st_sd, sr3c_vdem2_st_violsd, srg3c_vdem2_st_sd, srg3c_vdem2_st_violsd,se=c(cse_sr3c_vdem2_st_sd, cse_sr3c_vdem2_st_violsd, cse_srg3c_vdem2_st_sd, cse_srg3c_vdem2_st_violsd), column.labels = c("Model 1", "Model 2", "Model 3", "Model 4"),omit=c("years_no","yearf","factor"),  dep.var.labels.include=TRUE, dep.var.labels = c("SDM onset","Violent SDM onset","SDM onset","Violent SDM onset"), style = "ajps", notes.append = TRUE, notes.align = "l", model.numbers = F, notes = "Country-clustered errors in parentheses. Cubic term for time dependence of dependent variable and year- and region-fixed effects included but not reported.", order=c(paste0("^", vars_vdem2, "$"), ":", paste0("^", vars.order_main2, "$")), covariate.labels = c("Autonomy","Transition (0-5)","Autonomy x transition (0-5)","Transition proximity","Autonomy x transition proximity","Most powerful","Included","Group size","TEK irredentism","TEK SDMs","Petroleum % in settlement area","Previous SDMs in country","Normalized polity score","Distance border (log)","GDP p.c. (log)","Population (log)","Ongoing SD"))
###b) Wald tests
sr3c_vdem2_st_sd_jointsig <- linearHypothesis(sr3c_vdem2_st_sd, "autonomy + d5_tt_vdem2:autonomy = 0", test = "Chisq", vcov = cluster.vcov(sr3c_vdem2_st_sd, subset(df_global, sd_l1 == 0 & d5_tt_vdem2_st == 0)$cowcode))
sr3c_vdem2_st_sd_equality <- linearHypothesis(sr3c_vdem2_st_sd, "autonomy = autonomy + d5_tt_vdem2:autonomy", test = "Chisq", vcov = cluster.vcov(sr3c_vdem2_st_sd, subset(df_global, sd_l1 == 0 & d5_tt_vdem2_st == 0)$cowcode))
sr3c_vdem2_st_violsd_jointsig <- linearHypothesis(sr3c_vdem2_st_violsd, "autonomy + d5_tt_vdem2:autonomy = 0", test = "Chisq", vcov = cluster.vcov(sr3c_vdem2_st_violsd, subset(df_global, violsd_l1 == 0 & d5_tt_vdem2_st == 0)$cowcode))
sr3c_vdem2_st_violsd_equality <- linearHypothesis(sr3c_vdem2_st_violsd, "autonomy = autonomy + d5_tt_vdem2:autonomy", test = "Chisq", vcov = cluster.vcov(sr3c_vdem2_st_violsd, subset(df_global, violsd_l1 == 0 & d5_tt_vdem2_st == 0)$cowcode))
srg3c_vdem2_st_sd_jointsig <- linearHypothesis(srg3c_vdem2_st_sd, "autonomy + tt_vdem2_shock:autonomy = 0", test = "Chisq", vcov = cluster.vcov(srg3c_vdem2_st_sd, subset(df_global, sd_l1 == 0 & d5_tt_vdem2_st == 0)$cowcode))
srg3c_vdem2_st_sd_equality <- linearHypothesis(srg3c_vdem2_st_sd, "autonomy = autonomy + tt_vdem2_shock:autonomy", test = "Chisq", vcov = cluster.vcov(srg3c_vdem2_st_sd, subset(df_global, sd_l1 == 0 & d5_tt_vdem2_st == 0)$cowcode))
srg3c_vdem2_st_violsd_jointsig <- linearHypothesis(srg3c_vdem2_st_violsd, "autonomy + tt_vdem2_shock:autonomy = 0", test = "Chisq", vcov = cluster.vcov(srg3c_vdem2_st_violsd, subset(df_global, violsd_l1 == 0 & d5_tt_vdem2_st == 0)$cowcode))
srg3c_vdem2_st_violsd_equality <- linearHypothesis(srg3c_vdem2_st_violsd, "autonomy = autonomy + tt_vdem2_shock:autonomy", test = "Chisq", vcov = cluster.vcov(srg3c_vdem2_st_violsd, subset(df_global, violsd_l1 == 0 & d5_tt_vdem2_st == 0)$cowcode))

##########6.3 Excluding groups with highest power status##########
###a) models
sr10_vdem2_sd <- glm(as.formula(paste("sd_onset", paste(c("d5_tt_vdem2", "autonomy", "autonomy:d5_tt_vdem2","included","yearf","as.factor(region)","included","size","log(distance_border+0.00001)","tek_sup_irre_any","tek_sdm","oil_area_pct",country_controls, sdm_years), collapse = " + "), sep = " ~ ")), data=subset(df_global, state_control == 0 & sd_l1 == 0), family=binomial(link="logit"),na.action = na.exclude)
cse_sr10_vdem2_sd <- data.frame(cluster.se(sr10_vdem2_sd,as.factor(subset(df_global, state_control == 0 & sd_l1 == 0)$cowcode))[, 2])
sr10_vdem2_violsd <- glm(as.formula(paste("violsd_onset", paste(c("d5_tt_vdem2", "autonomy", "autonomy:d5_tt_vdem2","included","yearf","as.factor(region)","included","size","log(distance_border+0.00001)","tek_sup_irre_any","tek_sdm","oil_area_pct",country_controls, violsdm_years), collapse = " + "), sep = " ~ ")), data=subset(df_global, state_control == 0 & violsd_l1 == 0), family=binomial(link="logit"),na.action = na.exclude)
cse_sr10_vdem2_violsd <- data.frame(cluster.se(sr10_vdem2_violsd,as.factor(subset(df_global, state_control == 0 & violsd_l1 == 0)$cowcode))[, 2])
srg10_vdem2_sd <- glm(as.formula(paste("sd_onset", paste(c("tt_vdem2_shock*autonomy","included","yearf","as.factor(region)","included","size","log(distance_border+0.00001)","tek_sup_irre_any","tek_sdm","oil_area_pct",country_controls, sdm_years), collapse = " + "), sep = " ~ ")), data=subset(df_global, state_control == 0 & sd_l1 == 0), family=binomial(link="logit"),na.action = na.exclude)
cse_srg10_vdem2_sd <- data.frame(cluster.se(srg10_vdem2_sd,as.factor(subset(df_global, state_control == 0 & sd_l1 == 0)$cowcode))[, 2])
srg10_vdem2_violsd <- glm(as.formula(paste("violsd_onset", paste(c("tt_vdem2_shock*autonomy","included","yearf","as.factor(region)","included","size","log(distance_border+0.00001)","tek_sup_irre_any","tek_sdm","oil_area_pct",country_controls, violsdm_years), collapse = " + "), sep = " ~ ")), data=subset(df_global, state_control == 0 & violsd_l1 == 0), family=binomial(link="logit"),na.action = na.exclude)
cse_srg10_vdem2_violsd <- data.frame(cluster.se(srg10_vdem2_violsd,as.factor(subset(df_global, state_control == 0 & violsd_l1 == 0)$cowcode))[, 2])
stargazer(type="text",sr10_vdem2_sd, sr10_vdem2_violsd, srg10_vdem2_sd, srg10_vdem2_violsd,se=c(cse_sr10_vdem2_sd, cse_sr10_vdem2_violsd, cse_srg10_vdem2_sd, cse_srg10_vdem2_violsd), column.labels = c("Model 1", "Model 2", "Model 3", "Model 4"),omit=c("years_no","yearf","factor"),  dep.var.labels.include=TRUE, dep.var.labels = c("SDM onset","Violent SDM onset","SDM onset","Violent SDM onset"), style = "ajps", notes.append = TRUE, notes.align = "l", model.numbers = F, notes = "Country-clustered errors in parentheses. Cubic term for time dependence of dependent variable and year- and region-fixed effects included but not reported.", order=c(paste0("^", vars_vdem2, "$"), ":", paste0("^", vars.order_main2, "$")), covariate.labels = c("Autonomy","Transition (0-5)","Autonomy x transition (0-5)","Transition proximity","Autonomy x transition proximity","Included","Group size","TEK irredentism","TEK SDMs","Petroleum % in settlement area","Previous SDMs in country","Normalized polity score","Distance border (log)","GDP p.c. (log)","Population (log)","Ongoing SD"))
###b) Wald tests
sr10_vdem2_sd_jointsig <- linearHypothesis(sr10_vdem2_sd, "autonomy + d5_tt_vdem2:autonomy = 0", test = "Chisq", vcov = cluster.vcov(sr10_vdem2_sd, subset(df_global, state_control == 0 & sd_l1 == 0)$cowcode))
sr10_vdem2_sd_equality <- linearHypothesis(sr10_vdem2_sd, "autonomy = autonomy + d5_tt_vdem2:autonomy", test = "Chisq", vcov = cluster.vcov(sr10_vdem2_sd, subset(df_global, state_control == 0 & sd_l1 == 0)$cowcode))
sr10_vdem2_violsd_jointsig <- linearHypothesis(sr10_vdem2_violsd, "autonomy + d5_tt_vdem2:autonomy = 0", test = "Chisq", vcov = cluster.vcov(sr10_vdem2_violsd, subset(df_global, state_control == 0 & violsd_l1 == 0)$cowcode))
sr10_vdem2_violsd_equality <- linearHypothesis(sr10_vdem2_violsd, "autonomy = autonomy + d5_tt_vdem2:autonomy", test = "Chisq", vcov = cluster.vcov(sr10_vdem2_violsd, subset(df_global, state_control == 0 & violsd_l1 == 0)$cowcode))
srg10_vdem2_sd_jointsig <- linearHypothesis(srg10_vdem2_sd, "autonomy + tt_vdem2_shock:autonomy = 0", test = "Chisq", vcov = cluster.vcov(srg10_vdem2_sd, subset(df_global, state_control == 0 & sd_l1 == 0)$cowcode))
srg10_vdem2_sd_equality <- linearHypothesis(srg10_vdem2_sd, "autonomy = autonomy + tt_vdem2_shock:autonomy", test = "Chisq", vcov = cluster.vcov(srg10_vdem2_sd, subset(df_global, state_control == 0 & sd_l1 == 0)$cowcode))
srg10_vdem2_violsd_jointsig <- linearHypothesis(srg10_vdem2_violsd, "autonomy + tt_vdem2_shock:autonomy = 0", test = "Chisq", vcov = cluster.vcov(srg10_vdem2_violsd, subset(df_global, state_control == 0 & violsd_l1 == 0)$cowcode))
srg10_vdem2_violsd_equality <- linearHypothesis(srg10_vdem2_violsd, "autonomy = autonomy + tt_vdem2_shock:autonomy", test = "Chisq", vcov = cluster.vcov(srg10_vdem2_violsd, subset(df_global, state_control == 0 & violsd_l1 == 0)$cowcode))

##########6.4 Linear probability model##########
###a) models
sr3_vdem2_sd <- lm(as.formula(paste("sd_onset", paste(c("d5_tt_vdem2", "autonomy", "autonomy:d5_tt_vdem2","state_control","included","yearf","as.factor(region)",group_controls,country_controls, sdm_years), collapse = " + "), sep = " ~ ")), data=subset(df_global, sd_l1 == 0),na.action = na.exclude)
cse_sr3_vdem2_sd <- data.frame(cluster.se(sr3_vdem2_sd,as.factor(subset(df_global, sd_l1 == 0)$cowcode))[, 2])
sr3_vdem2_violsd <- lm(as.formula(paste("violsd_onset", paste(c("d5_tt_vdem2", "autonomy", "autonomy:d5_tt_vdem2","state_control","included","yearf","as.factor(region)",group_controls,country_controls, violsdm_years), collapse = " + "), sep = " ~ ")), data=subset(df_global, violsd_l1 == 0),na.action = na.exclude)
cse_sr3_vdem2_violsd <- data.frame(cluster.se(sr3_vdem2_violsd,as.factor(subset(df_global, violsd_l1 == 0)$cowcode))[, 2])
srg3_vdem2_sd <- lm(as.formula(paste("sd_onset", paste(c("tt_vdem2_shock*autonomy","state_control","included","yearf","as.factor(region)",group_controls,country_controls, sdm_years), collapse = " + "), sep = " ~ ")), data=subset(df_global, sd_l1 == 0),na.action = na.exclude)
cse_srg3_vdem2_sd <- data.frame(cluster.se(srg3_vdem2_sd,as.factor(subset(df_global, sd_l1 == 0)$cowcode))[, 2])
srg3_vdem2_violsd <- lm(as.formula(paste("violsd_onset", paste(c("tt_vdem2_shock*autonomy","state_control","included","yearf","as.factor(region)",group_controls,country_controls, violsdm_years), collapse = " + "), sep = " ~ ")), data=subset(df_global, violsd_l1 == 0),na.action = na.exclude)
cse_srg3_vdem2_violsd <- data.frame(cluster.se(srg3_vdem2_violsd,as.factor(subset(df_global, violsd_l1 == 0)$cowcode))[, 2])
stargazer(type="text",sr3_vdem2_sd, sr3_vdem2_violsd, srg3_vdem2_sd, srg3_vdem2_violsd,se=c(cse_sr3_vdem2_sd, cse_sr3_vdem2_violsd, cse_srg3_vdem2_sd, cse_srg3_vdem2_violsd), column.labels = c("Model 1", "Model 2", "Model 3", "Model 4"),omit=c("years_no","yearf","factor"),  dep.var.labels.include=TRUE, dep.var.labels = c("SDM onset","Violent SDM onset","SDM onset","Violent SDM onset"), style = "ajps", notes.append = TRUE, notes.align = "l", model.numbers = F, notes = "Country-clustered errors in parentheses. Cubic term for time dependence of dependent variable and year- and region-fixed effects included but not reported.", order=c(paste0("^", vars_vdem2, "$"), ":", paste0("^", vars.order_main2, "$")), covariate.labels = c("Autonomy","Transition (0-5)","Autonomy x transition (0-5)","Transition proximity","Autonomy x transition proximity","Most powerful","Included","Group size","TEK irredentism","TEK SDMs","Petroleum % in settlement area","Previous SDMs in country","Normalized polity score","Distance border (log)","GDP p.c. (log)","Population (log)","Ongoing SD"))
###b) Wald tests
sr3_vdem2_sd_jointsig <- linearHypothesis(sr3_vdem2_sd, "autonomy + d5_tt_vdem2:autonomy = 0", test = "Chisq", vcov = cluster.vcov(sr3_vdem2_sd, subset(df_global, sd_l1 == 0)$cowcode))
sr3_vdem2_sd_equality <- linearHypothesis(sr3_vdem2_sd, "autonomy = autonomy + d5_tt_vdem2:autonomy", test = "Chisq", vcov = cluster.vcov(sr3_vdem2_sd, subset(df_global, sd_l1 == 0)$cowcode))
sr3_vdem2_violsd_jointsig <- linearHypothesis(sr3_vdem2_violsd, "autonomy + d5_tt_vdem2:autonomy = 0", test = "Chisq", vcov = cluster.vcov(sr3_vdem2_violsd, subset(df_global, violsd_l1 == 0)$cowcode))
sr3_vdem2_violsd_equality <- linearHypothesis(sr3_vdem2_violsd, "autonomy = autonomy + d5_tt_vdem2:autonomy", test = "Chisq", vcov = cluster.vcov(sr3_vdem2_violsd, subset(df_global, violsd_l1 == 0)$cowcode))
srg3_vdem2_sd_jointsig <- linearHypothesis(srg3_vdem2_sd, "autonomy + tt_vdem2_shock:autonomy = 0", test = "Chisq", vcov = cluster.vcov(srg3_vdem2_sd, subset(df_global, sd_l1 == 0)$cowcode))
srg3_vdem2_sd_equality <- linearHypothesis(srg3_vdem2_sd, "autonomy = autonomy + tt_vdem2_shock:autonomy", test = "Chisq", vcov = cluster.vcov(srg3_vdem2_sd, subset(df_global, sd_l1 == 0)$cowcode))
srg3_vdem2_violsd_jointsig <- linearHypothesis(srg3_vdem2_violsd, "autonomy + tt_vdem2_shock:autonomy = 0", test = "Chisq", vcov = cluster.vcov(srg3_vdem2_violsd, subset(df_global, violsd_l1 == 0)$cowcode))
srg3_vdem2_violsd_equality <- linearHypothesis(srg3_vdem2_violsd, "autonomy = autonomy + tt_vdem2_shock:autonomy", test = "Chisq", vcov = cluster.vcov(srg3_vdem2_violsd, subset(df_global, violsd_l1 == 0)$cowcode))



